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ABSTRACT 

w' 

p^ I Simultaneous, precise measurements of the mass M and radius R of neutron stars can yield 

,^ '_ uniquely valuable information about the still uncertain properties of cold matter at several 

times the density of nuclear matter. One method that could be used to measure M and 
R is to analyze the energy-dependent waveforms of the X-ray flux oscillations seen during 



Ph! 



O 
H 
j/j ■ some thermonuclear bursts from some neutron stars. These oscillations are thought to be 

^ produced by X-ray emission from hot regions on the surface of the star that are rotating at or 

near the spin frequency of the star. Here we explore how well M and R could be determined 

^ I by analyzing energy-resolved X-ray data obtained using a future space mission having 2-30 

^^ I keV energy coverage and an effective area of 10 m^, such as the proposed LOFT or AXTAR 

^^ ' missions. We do this by generating energy-dependent synthetic observed waveforms for a 

CN ■ variety of neutron star and hot spot properties and then using a Bayesian approach and 

■^ ■ Markov chain Monte Carlo sampling methods to determine the joint posterior probability 

^f-^ . distributions of the parameters in our waveform model, given each synthetic waveform. We 

use the resulting posterior distributions to determine Bayesian confidence regions in the 

M-R plane by marginalizing the other parameters in our model. We explore how the sizes 

^ I and positions of these confidence regions depend on the inclinations of the hot spot and the 

observer, the background count rate, and deviations in the actual shape of the hot spot, 

radiation beaming pattern, and spectrum from those assumed in the waveform model. We 

also explore the effect on the confidence regions if the distance to the star or the inclination 

of the observer are known from other measurements, if a resonance scattering line is observed 

in the burst oscillation spectrum, or if the properties of the background are independently 

known. We assume that about 10^ counts are collected from the hot spot and that all sources 

of background contribute about 0.3 x 10^, 10^, or 9 x 10^ counts. 

We find that the uncertainties in the measured values of M and R depend strongly on 

the inclination of the hot spot relative to the spin axis. If the hot spot is within 10° of the 

rotation equator, both M and R can usually be determined with an uncertainty of about 

10%. If instead the spot is within 20° of the rotation pole, the uncertainties are so large 

that waveform measurements alone provide no useful constraints on M and R. Observation 

of an identifiable atomic line in the hot-spot emission always tightly constrains M/R; it 

can also tightly constrain M and R individually if the spot is within about 30° of the 



S 



rotation equator. These constraints can usually be achieved even if the burst oscillations 
vary with time and data from multiple bursts must be used to obtain 10® counts from 
the hot spot. Independent knowledge of the observer's inclination can greatly reduce the 
uncertainties, as can independent information about the background. Knowledge of the 
star's distance can also help, but not as much. Modest deviations of the actual spectrum 
from that assumed in the waveform model have little effect on the accuracies or uncertainties 
of M and R estimates. Large deviations of the actual shape of the hot spot or the radiation 
beaming pattern from those assumed can sometimes increase the uncertainties significantly, 
and can in some cases bias M and R estimates by moderate amounts. Our results show 
that if a sufficient number of burst oscillations produced by hot spots at high inclinations 
are observed using the next generation of X-ray timing satellites, one can constrain M and 
R tightly enough to discriminate strongly between competing models of cold, high-density 
matter. 

Subject headings: stars: neutron — stars: rotation — X-rays: bursts — X-rays: stars 



INTRODUCTION 



The properties of matter at extremely high densities are among the most important currently 
unresolved questions in physics and astronomy. Neutron stars contain large quantities of cold matter 
at densities that are otherwise inaccessible. Studies of these stars can therefore help determine the 
properties of such matter. In particular, simultaneous measurements of the mass M and radius R of 
neutron stars could provide tight constraints on the equation of state of ultradense matter (see, e .g. 
Lattimerll2nn7l : iLattimer fc PrakashI 120071 : IPead et al.lbood : lOzel fc Psalti^bood : iHebeler et aPboiol ). 



Following the discov ery of thermonucl ear X-ray bursts fro m neutron stars in the mid-1970s (jGrindlav et al 



19761 : iLewin et al.lll976l : for a review, see iLewin et al.lll993l ) and the discovery two decades later that 



some of these bursts pro duce X-ray f lux oscillations at or near the star's spin frequency (jStrohmaver et al. 



19961 : for a review, see IWattsll2012l ). several approaches have been proposed for using observations of 
X-ray bursts to determine M and R. 

The first attempts to constrain ne utron star properti e s using X-ray bursts used their apparent pea k 



luminosity and spectral temperature (jvan Paradiis 



1979 



Goldmanlll979l : lHoshi 



1981 



Marshall llMI)- 



With improved understanding of th e complexities of b urst emission, this approach was developed into a 
method for measuring M and R (see lLewin et al.lll993l ). This method relies on five assumptions: (1) that 
the radius of the photosphere at "touchdown" (defined as the moment after the photosphere has reached 
its maximum inferred radius when the temperature Tbb derived from fitting a Planck spectrum to the 
observed spectrum is highest) is the stellar radius; (2) that the emitting area at "touchdown" is the 
entire surface of the star; (3) that the flux seen at Earth at touchdown is the Eddington flux diluted by 
the distance to the burst source; (4) that the distance to the source is known; and (5) that the ratio of 
Tbb to the effective temperature T^s is known. 



This first method has been widely used to con strain M and R (see, e.g., Ivan Paradijsl 1 19821 : 
Paczyhski &: AndersonI 119861 : Ivan Paradijs et al.l Il990l ) . Tight constraints on M and R have recently 
been derived by us ing this method to analyze X-ray burst data obtained using the Rossi X-ray Timing 
Explorer (RXTE) (I6zelll2006l : l6zel et ahlliooal . boiol : lOiiver et al.lboiOal lblv The derived constraints on 



M and R are substantially smaller than the uncertainties in the observed values of the input parameters 
because most of th e observed values are disc arded, since they wo uld produce values of M and R that are 
complex numbers (jOzel et al.1 120101 : see also ISteiner et alj l20ld ) . When this inconsistency is addressed 
and systematic errors are fully inc luded, the uncertainties in M and R are likely to be substantially 
larger (see, e.g.. ISteiner et al.ll2010l ). 



A second approach to measuring M and R using X-ray burs t observations is to fit det a iled spectral 



models to high-precisio n measurements of X-ray burst spectra (jMaiczvna &: Madeil l2005l : iMiller et al 



201 ll : iMiller et al.ll2013l ) . This approach assumes that the spectrum predicted by the atmospheric model 
is an accurate physical description of the observed spectrum. If it is, fitting a grid of spectral models 
to the observed X-ray spectra yields the radiation temper ature Tm in the l o cally comoving frame a t 



the stellar surface as well as M, R, and the composition (JMiller et al.ll201ll : ISuleimanov et al.ll2012l ). 
This approach can only be used if high-precision spectral data and accurate, high-precision spectral 
models are both available. At presen t, adequate data are available only for the 4U 1820—30 superburst 



(JMiller et al. 



2011 



Miller et al.ll201.'j ). 



A third approach for determining M and R using X- ray burst observations is to fit models of burst 
oscillati on waveforms to obser vatioi is of these wave f orms (jStrohmaver et al.lll997l : IMiller &: Lamblll998l : 
see also [Weinberg et al.ll200ll ). As IMiller &: LambI ( 19981 ) demonstrated, fitting energy-resolved wave- 
form data can improve the constraints. Burst oscillations are thought to be produced by X-ray emission 
from a region on the surface of the sta r that is hotter than th e rest of the surf ace and is ro tating at 
or near the spin frequency of the star (jStrohmayer et al.lll996l : for a review, see IWattsI l2012l ). Such a 
hotter region could be present either because a part of the stellar surface has retained more heat than 
other parts after thermonuclear burning has occurred or because a disturbance in the outer layers of 
the star, such as a surface normal mode, has made a localized region hotter. The observed amplitude 
of the waveform constrains the inclination of the hot spot relative to the spin axis, the observer's incli- 
nation, and the compactness (M/R) of the star, the last primarily via general relativistic light-bending 
effects. The observed asymmetry and harmonic content of the waveform constrain the component of 
the velocity of the emitting region in the observer's direction, primarily via special relativistic Doppler 
boosts and aberration. Because the rotation frequency of the emitting region is accurately known from 
the oscillation frequency, knowing the surface velocity constrains the stellar radius. 



Nath et al.l (|2002l ) explored the constraints that could be derived on the compactness of 4U 1636—536 



and 4U 1728—34 by analyzing RXTE observations of the bolometric fiux oscillations that occur during 
the rise of X-ray bursts from these two neutron stars. They modeled the bolometric waveforms during 
the burst rise using small hot spots that expand linearly with time, neglected the relativistic Doppler 
shifts and aberration produced by the rotational motion of the hot spots, and assumed that the back- 
ground was known. They found that they could not determine whether the oscillations are produced by 
a single spot or two antipodal spots, or constrain the hot spot and viewing geometry assuming either 
of these alternatives, using bolometric RXTE waveform data. The primary reason is that the changes 
in the waveform produced by changes in different parameters of the model are very similar. As a result 
of this degeneracy, they were unable to constrain the stellar compactness for models with a single hot 
spot, even if they assumed that the spot and the observer were known from other information to be in 
the rotation equator. They were able to obtain interesting upper bounds on the compactness for models 
with two antipodal hot spots and argued that it would be possible to simultaneously constrain the stellar 
compactness and the hot spot and viewing geometries with a count rate ~ 10-20 times higher than the 
RXTE rate. However, the 290 Hz subharmonic in the 4U 1636—536 waveform that appeared significant 



in one burst (|MilleiJ Il999l : IStrohmaveiJ 120011 ) , and suggested consideration of two antipod al spots, was 
subse quently found not to be significant in a more detailed analysis of additional data ([Strolimayer 



2001 



Bhattacharyya et alJ (120051 ) investigated the constraints on the stellar compactness, hot spot prop- 
erties, and system geometry that could be obtained by fitting model pulse profiles to an average energy- 
resolved oscillation profile produced by folding and stacking RXTE observations of 22 burst oscillations 
in three groups of bursts fro m XTE J1814— 338. Sta ble oscillations occur throughout the bursts pro- 
duced by this neutron star (IStrohniaver et al.1 12003| ) . In generating their energy-resolved oscillation 
profile models, iBhattacharvva et alj assumed a single hot spot and an oscillation-phase-independent 
background chosen so that the model produces the observed total number of counts in each energy 
channel, when summed over oscillation phase. They included the relativistic Doppler shifts, aberration, 
and frame dragging produced by the rotation of t he star, as well as gr avitational redshift and light 
bending effects. In order to include frame dragging, IBhattacharvva et al.l had to construct general rela- 
tivistic stellar models numerically and therefore considered just two illustrative equations of state. They 
were able to achieve acceptable fits to the average oscillation profile using these equations of state, but 
obtained only weak constraints on the stellar compactness, hot spot properties, and system geometry. 



Strohmaveii (|2004l ) explored the constraints on M and R that might be achieved by analyzing burst 



oscillation observations with much higher count rates than were achieved using RXTE. H e did this by 



analyzing a synthetic bolometric waveform similar to a waveform seen in 4U 1636—536. IStrohmavei 



NathetaL 



genera ted his synthetic and model waveforms using a code that is similar to the one used by 
(I2OO2I ) but included the relativistic Doppler shifts and aberration produced by the rotation al motion of 
the hot spots. The synthetic observed waveform did not include a background of any kind. IStrohmayer 
then assumed that the absence of any background is known to the observer by means other than fitting 
the waveform]^ He also assumed that the spot and the observer are both in the rotation equator, which 
is the most favorable possible geometry, and made the strong assumptions that the hot spot size and 
location, the radiation spectrum and beaming pattern, and the observer's inclination are known inde- 
pendently of the waveform, so that only M and R have to be estimated using the waveform data. These 
assumptions eliminate the strong degeneracies between M, i?, and the other model par ameters that 
greatly increase the uncertainties in more realistic situations. Based on these assumptions. IStrohmayer 
concluded that an X-ray timing mission with a collecting area ~ 10 times larger than RXTE would 
be able to determine M and R to within a few percent, thereby placing interesting constraints on the 
equation of state of neutron star matter. 



Muno et al.l (|2003l ) examined averages of the energy-resolved flux oscillations observed from several 



X-ray burst sources using RXTE. They found that folded oscillations observed in their higher energy 
bands arrived later than those in their lower energy bands, although the energy dependence varied 
significantly with epoch and source. Their analysis of folded and averaged 4U 1636—536 oscillation 



^We note that attempting to remove the background by subtracting the pre-burst count rate assumes that the pre- 
burst background, which is presumably produced at least in part by accretion onto the star, persists unchanged throughout 
the burst. This introd uces possible systematic errors, because luminous bursts are expected to alter the accretion flow 
significantly (see, e.g.. iMiller fc Lamblll996l : IWorpel et al.ll2013h . Even if the pre-burst background persists unchanged 
throughout the burst, subtracting it from the count rate during the burst, rather than including the background as a 
component of the model, incorrectly neglects the fluctuation in the number of counts produced by the background and the 
corresponding uncertainties in the model parameters. This conceptual error is built into the current version of the XSPEC 
analysis package, so all data analysis done using XSPEC makes this error. 



profiles showed the clearest evidence for suc h a trend, which wou ld be inconsistent with a simple rotating 
spot model of burst oscillations. Recently, lArtigue et alj (|2013l ) have analyzed the same 4U 1636—536 
data in much more detail and find that the data are entirely consistent with a simple rotating spot 
model, although the parameters of such models are poorly constrained by these data. 

New mission concepts are now being proposed that would provide much larger count rates than pre- 
vious missions. Plans for the proposed Large Observatory for X-ray T iming (LOFT) include a detector 
with an effective area ~ lOm^ at 8 keV and 2-30 keV energy coverage ( Feroci et al.ll201Cll : lMignani et al 



2012l : lDel Monte et al.ll2012l ). more than an order of magnitude larger than the area of the RXTE Pro- 
portional Counter Array (PCA). The Advanced X-ray Timing Array [AXTAR) concept includes a 



detector with an effe ctive area greater than 3 m^ and 2-50 keV energy coverage ( Chakrabartv et al 



20081 : iRav et al.ll201lh . One of the main motivations for these missions is to determine the properties of 
neutron stars with much higher precision than has been possible previously, by analyzing high-precision 
observations of the waveforms of burst oscillations. Although we focus here on exploring the constraints 
on M and R that can be obtained by analyzing the waveforms of X-ray burst oscillations, our results 
are equally relevant for measuring M and R using the waveforms produced b y X-ray emissioii frorn . 
the h eated polar caps of isolated rotation-powered millis econd pulsars (see, e. g. iBogdanov et al.l 120071 . 
20081 ). w hich is the goal o f the proposed NICER mission (JGendreau et al.ll2012l ). Using methods similar 



to ours, iBogdanovl (|2013l ) has recently derived a lower limit of 11.1 km on the radius of the isolated 
pulsar PSR J0437— 4715, assuming its mass is 1.76 Mq and that systematic errors can be neglected. 

In this paper we explore the constraints on M and R that could be derived by analyzing energy- 
resolved burst oscillation waveforms obtained using a future, satellite-borne detector with 2-30 keV 
energy coverage and an effective area 10 to 20 times larger than the RXTE PCA. We do this by first 
generating energy-dependent synthetic observed waveforms for a variety of neutron star and hot spot 
properties. We then use a Bayesian approach and Markov chain Monte Carlo (MCMC) sampling meth- 
ods to determine the constraints on M and R that can be obtained by analyzing these synthetic observed 
waveforms. Specifically, we determine the joint posterior probability distribution of the parameters in 
our waveform model, given the synthetic waveform of interest, use this distribution to determine the 
joint posterior distribution of M and R by marginalizing the other parameters in the waveform model, 
and then use the joint distribution of M and R to determine the most probable values of M and R and 
Bayesian confidence regions in the M-R plane. 

Our analysis applies to any waveforms produced by emission from a single region of hotter gas that 
is rotating about the star. Such a region could be present either because a part of the stellar surface has 
retained more heat than other parts after thermonuclear burning has occurred or because a disturbance 
in the outer layers of the star, such as a surface normal mode, has made a localized region hotter (see 
Wattsll2012l ). Oscillations with the amplitudes ~ 10% that are required to obtain significan t constraints 
on M and R are very probably produced by a single hotter region (see lLamb et al.ll2009al ). 



We explore how the sizes and positions of these confidence regions depend on the inclinations of the 
hot spot and the observer and the background count rate. We also explore the effect on the confidence 
regions if the distance to the star or the inclination of the observer are known from other measurements, 
if a resonance scattering line is observed in the burst oscillation spectrum, or if the properties of the 
background are independently known. Finally, we explore the effects of deviations in the actual shape 
of the hot spot, radiation beaming pattern, and spectrum from those assumed in the fitted model. 
We assume that about 10^ counts are collected from the hot spot and that all sources of background 



contribute about 0.3 x 10^, 10^, or 9 x 10^ counts. Our treatment of the background is very conservative, 
in the sense that we usuahy make no assumptions about the magnitude or spectrum of the background. 
We do not even assume that the background is constant, only that it does not vary at frequencies 
commensurate with the burst oscihation frequency. 

We find that the uncertainties in the measured values of M and R depend strongly on the inclination 
of the hot spot relative to the stellar rotation axis. If the hot spot is within 10° of the rotation equator, 
both M and R can usually be determined with an uncertainty of about 10%. If instead the spot is 
within 20° of the rotation pole, the uncertainties are so large that waveform measurements alone provide 
no useful constraints on M and R. The uncertainties in M and R are affected little by background 
count rates less than or comparable to the count rate from the hot spot, but become significantly larger 
for higher background count rates. The precisions we report here can usually be achieved even if the 
burst oscillations vary with time and data from multiple bursts must be combined to obtain 10^ counts 
from the hot spot. 

Observation of an identifiable atomic line in the hot-spot emission always tightly constrains M/R; 
it can also tightly constrain M and R individually, if the spot is within about 30° of the rotation 
equator. Independent knowledge of the observer's inclination can greatly reduce the uncertainties, as 
can independent information about the background. Knowledge of the star's distance can also help, 
but not as much. 

Modest deviations of the actual spectrum from that assumed in the fitted model have little effect 
on the accuracy or uncertainty of M and R estimates. Large deviations of the actual radiation beaming 
pattern from the pattern assumed in the waveform model can increase the uncertainties of M and R 
measurements substantially. In some cases, but not always, large deviations of the actual shape of 
the hot spot from the circular shape assumed in the waveform model can increase the uncertainties of 
M and R estimates and bias them by moderate amounts. The physical conditions that produce tight 
constraints on M and R (relatively small spots far from the rotation pole) are the conditions in which 
the shape of the spot is unimportant. 

Our results show that if a sufficient number of burst oscillations produced by hot spots at high 
inclinations are observed using the next generation of X-ray timing satellites, one can constrain M and 
R tightly enough to discriminate strongly between competing models of cold, high-density matter. 

The remainder of this paper is organized as follows. In Section [21 we outline our approach to 
generating synthetic observed waveforms and producing the model waveforms that we fit to the synthetic 
waveform data. In Section [3l we describe the computational methods we used to produce synthetic 
waveform data and construct model waveforms, and the MCMC computational methods we used to 
determine the posterior probability distributions of the parameters in the model, given a synthetic 
waveform. In Section [H we describe our results and in Section [5l we summarize our conclusions. 
In Appendix [Aj we summarize the suite of test problems and solutions that we use to validate our 
waveform and MCMC codes. In Appendix [HI we discuss the constraints on system parameters that can 
be obtained by jointly fitting many segments of waveform data, from a single burst or from multiple 
bursts. 



2. APPROACH 

2.1. Generation of "observed" and model waveforms 

In the present work we fit energy-resolved model waveforms to energy-resolved synthetic "observed" 
waveforms, using a Bayesian approach. For both waveforms, we assume that the oscillation is produced 
by emission from a single, uniformly emitting, hotter region on the stellar surface. Such a region could 
be present either because a part of the stellar surface has retained more heat than other parts after 
thermonuclear burning occ urred or bec ause a disturbance in the outer layers of the star has made a 



localized region hotter (see IWattsI I2OI2I ) 



We consider synthetic observed waveforms generated by circular and elongated hot spots (see 
Section l4.1.2p . but to reduce the number of fitted parameters our model waveforms assume the hot 
spot is circular (see Section [4. l.ip . For simplicity, we usually assume the hotter region emits radiation 
with a Planck spectrum and 100% efficiency but with the beaming pattern appropriate for an electron 
scattering atmosphere, assumptions that are mutually inconsistent. This inconsistency can be avoided 
when fitting real data by including an appropriate color factor in the emission model used to construct 
fitted waveforms. We also assume that the radiation propagating from the emitting area on the stellar 
surface reaches the observer without interacting with any ambient or intervening matter, but we do 
include background counts to illustrate the effects of possible background emission from the stellar 
surface, an accretion disk, and other sources in the field, as well as the instrumental background. 

In generating synthetic and model waveforms, we assume that the hotter region has a constant size 
and shape, is located at a fixed stellar rotational latitude, and rotates at a constant frequency. These 
assumptions are not as restrictive as they might at first appear. The reason is that constraints on M 
and R similar to those obtained for waveforms that satisfy these assumptions can usually be obtained 
by appropriately analyzing data from a single burst or from multiple bursts, provided the data contains 
the same number of counts as assumed in our analysis. In particular, constraints on M and R similar to 
those found here can usually be achieved even if the oscillation frequency and other physical parameters 
vary during the burst or from burst to burst (see Section 14.2.71 and Appendix [B]) . 

In constructing the synthetic waveforms, we include a constant background component. This 
component is a catch-all for all counts not produced by radiation from the hot spot. These counts 
could be produced by emission from unassociated sources in the field, the accretion disk, the non-spot 
portion of the star, instrumental backgrounds, or any combination of these. For simplicity, we model 
this background by adding emission from the entire stellar surface with the beaming pattern expected 
for an electron scattering atmosphere and a spectrum having the shape of a Planck spectrum with a 
temperature lower than the temperature of the hot spot. We normalize the background spectrum to 
achieve the desired number of background counts. The number of counts contributed by this emission 
is important, but not their detailed properties. 

In fitting model waveforms to the synthetic waveform, we treat the background component in 
the model waveform very conservatively, in the sense that we usually make no assumptions about its 
magnitude or spectrum. We do not even assume that the background is constant, only that it does 
not vary at frequencies commensurate with the hot spot rotation frequency. Any prior knowledge of 
the properties of the background can be used to restrict the background model and usually tighten the 
constraints on the values of the system parameters that can be derived from waveform observations. 



We compute the time- and energy-resolved waveforms that would b e seen by a distant ob server 



using the Schwarzschild plus Doppler (S-|-D) approximation introduced bv lMiller &: LambI (| 19981 ) . The 
S-|-D approximation treats exactly all special relativistic effects (such as relativistic Doppler boosts and 
aberration) produced by the rotational motion of the hot spot, but treats the star as spherical and uses 
the Schwarzschild spacetime to compute the general relativistic redshift, trace the propagation of light 
from the stellar surface to the observer, and calculate light travel-time effects. The S-l-D approximation 
does not include the effects of stellar oblateness or frame dragging. However, for the stars considered 
here, and indeed for any stars that do not both rotate rapidly and have very low compactness, the effects 
of stellar oblateness and f rame dragging are m inimal and are negligible compared to the uncertainties 



in the X-ray emission (see lCadeau et al.ll2007l ). 



The S-|-D approximation allows us to exploit the spherical symmetry of the spacetime to reduce 
substantially the number of integrations needed to trace rays from the emitting area to the observer. 
As a result, this code is much faster than codes that take into account the oblateness of the star, 
approximate the star's exterior spacetime using the Kerr spacetime, or use numerically computed stars 
and exterior spacetimes. This speed is essential, because finding the high-probability region of the 
parameter space of the waveform model and then determining the posterior probability density in this 
region with some precision requires accurate computation of the likelihood at a large number of points 
in a high-dimensional parameter space, and therefore requires many ray tracings. 

The waveform code we use is based on th e code used and val idated in the analysis of accretion- 



powered millisecond X-ray pulsar waveforms by iLamb et al.l (l2009al )bl) . We report the results of further 
code validation tests in Appendix lAl 

To determine the constraints on M and R that can be derived from burst oscillation waveforms, 
we compute the posterior probability distribution of all the parameters in the waveform model, for each 
of a variety of synthetic observed waveform data sets, using a standard Bayesian approach. We use an 
MCMC algorithm to sample the parameters of the waveform model and compute the likelihood of each 
set of parameters, given the synthetic waveform being considered. 

We construct each synthetic observed waveform data set by Poisson sampling the counts in each 
phase and energy bin of each synthetic observed waveform that we computed using the S-l-D approxi- 
mation, to mimic the statistical fluctuations that would be present in actual data. 

We generate the model waveforms that we compare with the synthetic waveform data using the 
same code that we use to generate the synthetic waveforms. The model waveforms include a possible 
background component with an arbitrary energy spectrum. 



2.2. Burst rise or burst tail? 

Previous interest in using X-ray burst oscillations to determine the masses and radii of neutron 



stars has focused largel y on using oscillations observed during the rise of bursts (see, e.g., iNath et al. 



2002 



Strohmayerj l2004l ) . primarily because of the large fractional modulation observed during the rise 
of some bursts and the evidence that emission early in the burst rise comes from a small, hotter 
region on the stellar surface. Although the oscillations observed during burst tails usually have smaller 
fractional amplitudes, tail oscillations usually last much longer than the oscillations during burst rises. 
We therefore consider the merits of using oscillations observed during burst tails as well as burst rises. 



Key factors to consider in evaluating the relative merits of using burst rise oscillations and burst tail 
oscillations include (a) the total number of oscillating counts, (b) the total number of counts collected 
during the observation, (c) the information contained in these counts, including any information encoded 
in the constant X-ray flux from the star, (d) the variation of the oscillation waveform with time during 
the measurement interval, and (e) the complexity of the emission pattern. 

We show here that for a next-generation timing instrument with a collecting area 20 times that 
of the RXTE PCA, the uncertainties in M and R estimates obtained by analyzing oscillations during 
the tails of bursts are likely to be as small or smaller than the uncertainties derived by analyzing the 
oscillations observed during the rise of bursts. Hence burst tail oscillations are likely to provide results 
that are equally good or even superior to the results obtained from burst rise oscillations. 



2.2.1. Statistical uncertainties in parameter estimates 

Consider first the uncertainties in M and R estimates produced by the fluctuations in the observed 
waveform caused by photon counting noise. We expect the fractional uncertainties in M and R estimates 
to decrease as the ratio TZ of the number of modulated counts to the fluctuation in the number of counts 
increases. For the burst oscillation waveforms, we use as a figure of merit the quantity 



n = N^sc/VKa = 1.4 /rms iVtot , (l) 

where A'osc is the number of counts in the oscillating component of the waveform, defined as the integral 
of the semi-amplitude of the oscillating count rate over the duration of the data segment; iVtot is the 
integral of the total count rate, including the background count rate, over the duration of the data 
segment; and /rms is the fractional rms amplitude of the oscillation during the data segment. Our 
Bayesian analysis (see Section 14.2. ip shows that the fractional sizes of the confidence regions in the 
M-R plane decrease with increasing 7^, approximately as 7^~^. The value of TZ is therefore a useful 
figure of merit when comparing different data sets. 

We illustrate the implications of this result for the relative merits of analyzing oscillations during 
burst rises and tails by scaling from the count rates observed during the rise and tail of a well-observed 
X-ray burst from 4U 1636—536 (oscillations were detected using RXT E in the rise of this burst but not 



in the tail). We f irst summarize the observed RXTE PCA count rates (jStrohmaver et al.lll998l : see also 



Nath et al.ll2002l ) and then use these count rates to estimate the count rates that might be observed 



using a next-generation large-area X-ray timing instrument like those mentioned in Section [TJ 

Oscillations during the burst rise. During the first 1/16 s of the burst, the average total count 
rate (including a background assumed equal to the ~2,000 counts s~^ pre-burst background) was 
~4, 500 counts s~^, while the average oscillation semi-amplitude was ~2,000 counts s~^. Thus, during 
the first 1/16 s of the burst a total of ~ 280 counts was collected, including ~ 125 modulated counts. 
Our postulated next-generation X-ray timing instrument would therefore collect ~ 5, 600 total counts, 
including ~ 2, 500 modulated counts, yielding an 7^- value ~ 33 when data from this part of the burst 
rise is used. 

During the first 1/4 s of the burst, which includes most of the burst rise, the average total count 
rate (including the assumed background) was ~8, 000 counts s~^, while the average oscillation semi- 
amplitude was ~ 2, 000 counts s~^ . Thus, during the first 1/4 s of the burst a total of ~ 2, 000 counts was 
collected, including ~ 500 modulated counts. Our postulated next-generation X-ray timing instrument 
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would therefore collect ~40, 000 total counts, including ~ 10, 000 modulated counts, yielding an 7^- value 
~50 when all the data from the burst rise is used. 

These results indicate that the 1/4 s data set would provide constraints similar to or even better 
than the 1/16 s data set, even though the average fractional amplitude of the oscillation was much 
smaller during the longer interval, because the total number of counts during the longer interval was 
much greater. 

Oscillations during the burst tail. During the last 5 s of the burst, the RXTE PCA collected a total 



of ~80, 000 counts. According to lGallowav et al.l (|2008l ). the mean fractional rms amplitude of a typical 



burst oscillation in the RXTE sample is ~ 10%, implying a semi-amplitude ~ 14% for a sinusoidal 
waveform. Our postulated next-generation X-ray timing instrument would therefore collect ~1.6 x 10^ 
total counts during the last 5 s of a similar burst, including ~200, 000 modulated counts if an oscillation 
with an rms amplitude of ~ 10% is present, yielding an 7^-value ~ 160. 

TZ scales as the square root of the number of counts, so combining data from the early part of 25 
bursts like this example could yield an 7^-value ~ 250, while combining data from the tails could yield 
an 7^- value ~800. 

These results suggest that analyses of data from X-ray burst tails may provide more precise con- 
straints than analyses of data from burst rises. Even if the tail oscillations last somewhat less than 5 s 
or the rms amplitude is somewhat less than ~ 10%, this example suggests that analyses of burst tail 
data may provide results comparable to those obtained by analyzing data from burst rises. 

In these estimates we have assumed that the unmodulated (constant in time) component of the 
X-ray flux time series contains no information about M and i?, i.e., that only the modulated component 
of the waveform contains such information. If instead the unmodulated component contains some infor- 
mation about the system (such as the compactness of the star, which tends to increase the unmodulated 
count rate due to increased gravitational lensing) and one can extract this information, then using the 
data collected during the burst tail could be more useful than suggested by the preceding estimates. 

To see this, suppose first that the X-ray waveform is constant in time (i.e., unmodulated) and 
that the magnitude of this constant count rate is important. Then the fractional uncertainties in 
estimates of M and R will depend only on the total number of counts. Now suppose that in addition 
to this same unmodulated component there is an infinitesimal modulated component. Clearly, this 
modulated component will not add much extra information about the system. This case may, in fact, 
apply to observations of the tails of some bursts. In this case, considering only the oscillations would 
underestimate our ability to constrain the model parameters using data from the burst tail compared 
to using data from the burst rise. 

These considerations suggest that using oscillations observed during the burst tail to estimate M 
and R is likely to be at least as effective as using oscillations observed during the burst rise, other things 
being equal. 



2.2.2. Other considerations 

In addition to the figure of merit TZ of the data set, other important considerations for assessing 
the relative merits of burst rise and burst tail oscillations include the time-dependence of the waveform 
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and the size and possibly complicated shape of the emitting region that produces the oscillations. Here 
we consider these two factors. 

Time- dependence of the waveform. The frequencies and amplitudes of tail oscillations usually 
change relatively slowly (see IWattsI l2012l ). In contrast, the oscillations seen during the rise of bursts 
often show large and rapid changes in amp litude and of ten in frequency, including substantial deviations 
from the stellar spin frequency (again see IWatta |2012| ). Hence, waveform models and procedures for 
fitting oscillations during the onset of bursts must be able to handle much more rapid and greater 
changes in the size, inclination, and longitude of the emitting region than are encountered in fitting tail 
oscillations. 

We have investigated the problem of analyzing burst oscillation waveforms that change with time 
(see Section [4.2.71 and Appendix iB]) . If the oscillation frequency is changing but the change can be 
modeled accurately enough to maintain the correct oscillation phase when folding successive periods of 
the oscillation, then the constraints on M and R that can be obtained by analyzing the resulting folded 
waveform will be nearly the same as those that could be obtained by analyzing a similar waveform 
with a fixed oscillation frequency and the same number of counts. If the oscillation frequency varies 
too rapidly or irregularly during the burst rise or tail to be described accurately by a simple frequency 
model, the full burst oscillation data set can be divided into smaller time segments and analyzed using 
standard Bayesian techniques. This approach can also be used if other physical properties of the system, 
such as the size and inclination of the emitting region, vary significantly. The computational burden of 
this kind of analysis increases only linearly with the number of segments. Consequently, variations of 
the burst oscillation waveform on timescales shorter than the burst rise or burst tail (but substantially 
longer than the burst oscillation period) do not appear to pose an insurmountable analysis problem. 

Complexity of the e mission pattern. Pr evious modeling of the waveforms of accretion-powered 
millisecond X-ray pulsars ( Lamb et al.l l2009al Jbl) showed that the size and shape of the emitting region 
have only a weak effect on the amplitude and shape of the waveform, unless the emitting region is very 
large (angular radius > 50° for inclinations > 60°). The hotter emitting regions that produ ce oscillation s 
early in the rise of a burst are thought to cover a small fraction of the stellar surface (jWattd l2012l ) . 
In this case the detailed shape of the emitting region does not significantly affect the properties of the 
waveform. Oscillations obs erved late in the rise of a burst or in burst tails are thought to be produced 
by larger emitting regions (jWattsll2012l ). but these regions may still be small enough that their detailed 
shape does not significantly affect the properties of the waveform. Our fits of waveform models that 
assume a circular emitting region to synthetic observed waveforms produced by large regions elongated 
in the east-west and north-south directions suggest that if the emitting region is large and distorted, 
this will increase the uncertainties in estimates of M and R, although not necessarily by a very large 
amount, but will not bias these estimates significantly (see Section l4.2.6p . Our experience in fitting 
waveforms suggests that results similar to those we report here are likely to be achievable, even if there 
are moderate temperature variations across the hot spot. 
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3. COMPUTATIONAL METHODS 

In this section we describe the computational methods we use to determine the accuracy and 
precision with which M and R could be determined using the energy- and time-resolved burst oscillation 
waveform data that could be obtained by future large-area X-ray timing missions. We apply standard 
Bayesian inference methods to compute the best-fit values and confidence intervals for the parameters 
M and R in our waveform model, for a variety of synthetic observed waveforms of interest. We first 
explain how we compute the oscillation waveforms that we use. We then discuss our Bayesian analysis 
and sampling methods. Finally, we describe how we integrate over uninteresting parameters to determine 
the joint posterior distribution of M and R and determine their best-fit values and confidence regions. 



3.1. Waveform Computation 

The burst oscillation waveform model that we use here has eight parameters. The neutron star is 
described by its gravitational mass M and circumferential radius R. We represent the emitting area 
as a circular, uniformly radiating hot spot on the stellar surface, with an angular radius A^gpot. We 
assume that the center of the spot is inclined at an angle ^spot relative to the star's spin axis and that 
the spot rotates uniformly around the spin axis of the star with a frequency f rot • We also assume that 
the spot emits radiation that has a blackbody spectrum with a temperature T^o when measured in an 
inertial frame at the stellar surface that is momentarily comoving with the surface. We assume that 
the values of these parameters are constant in time; the waveform is then perfectly periodic. Finally, 
we assume that the observer views the system from a distance d, at an inclination ^o^s relative to the 
star's spin axis. 

We use as our global coordinate system Schwarzschild coordinates (r, 6*, (/9, t) centered on the star, 
with = aligned with the spin axis and c/? = in the plane containing the spin axis and the observer. 
We choose the zero of the Schwarzschild time coordinate t so that a light pulse that propagates radially 
from a point on the stellar surface immediately below the observer (i.e., dX 9 = Oohs and y? = 0) arrives 
at the observer at t = 0. 

We use our waveform code to compute the phase- and energy-resolved photon number flux that 
would be seen by a distant observer during a single rotation of the hot spot. This code is based on the 



code w e used previously to compute the waveforms of accreting millisecond X-ray pulsars (jLamb et al. 
2009al )bl) and uses the Schwarzschild plus Doppler (S-l-D) approximation (JMiller fc Lamblll998l ) (see 



Section [2. ip . In this approximation, the exterior spacetime is specified completely by the stellar com- 
pactness, GM/Rc^ (hereafter, for conciseness we write M/R for GM/Rc^). The computational speed 
made possible by this approximation makes it practical to compute the required likelihood distributions 
of the model parameters on a large, modern computer cluster. 

It is convenient to specify the radiation from a point on the stellar surface in a local inertial frame 
located at the surface and momentarily comoving with it. We assume that the beaming pattern of 
the radiation emitted from the stellar surface is axisymmetric about the normal to the surface, when 
measured in the comoving frame. The specific intensity I'q in the comoving frame at the stellar surface 
can then be expressed as /g(£'Q,a'), where E'q and a' are the photon energy and the angle between the 
photon direction and the normal to the surface measured in the comoving frame at the stellar surface. 
We assume that Io(£'o,Q;') can be written as the product of a beaming function g{a') and a spectral 



- 13- 

function /{E'q), i.e., 

l',{Eia')=g{a')f{E',). (2) 

In the present work, we usually consider the waveforms produced by the beaming pattern expected 
for emission from an electron scattering atmosphere, but we sometimes consider isotropic beaming. For 
the former beaming function, we use the quadratic expression 

^'(a') = a + 6cosa' + ccos a , (3) 

with 

a = 0.42822, b = 0.92236, and c = -0.085751 . (4) 

This expression is a least-squares fit to the beaming function for emission from a uniform, semi-infinite, 
Thomson scattering atmospher e, taking polarization into account, and agrees with the actual beaming 



function (|Chandrasekhailll96d . Table XXIV) to better than 1% for 0.02 < cosq' < 1. 



In the present work, we usually set /'{Eq) equal to the Planck function B{Eq,T'), where T' is 
the radiation color temperature at the stellar surface, as measured in the comoving frame. We expect 
the actual spectral function to have a shape similar t o but not exactly the sa me as a Planck spectrum 



and an efficiency ~ 20%, rather than 100% (see, e.g.. ISuleimanov et al.ll2012l ). Using a color factor to 



produce the appropriate lower efficiency would be important in analyzing real data. 

In the few cases where we consider the presence in the spectrum of an atomic scattering line, we 
model the line by multiplying the continuum spectrum by the transmission factor exp[— r(£' — Ec)], 
where Ec is the centroid energy of the line and t(x) is a Gaussian profile having a maximum value r(0) 
and a specified FWHM. 

The specific intensity /q(£'q,q') measured in the comoving frame at the stellar surface can be 
converted to the specific intensity Io{Eq, a) measured in the static frame at the stellar surface using the 
invariance of I{E)/E'^ under a Lorentz boost. The result is 

/o(i^o,a) = /^(i?^,a')7"'[l " (v/c) • k]-^ , (5) 

where Eq is the photon energy, a is the angle between the unit vector k in the direction of the light ray 
and the normal to the surface, and v is the linear velocity of the stellar surface at the point of emission, 
all measured in the static frame, and 7 = [1 — (v/c)^]"^'^. The quantities £^0 f^d a are related to the 
corresponding quantities Eq and a' measured in the comoving frame by 

Eq = 6E'o (6) 

and 

cos a' = 6 cos a , (7) 

where 

S = ^ (8) 

7[l-(v/c)-k] 

is the Doppler factor. Here v • k = ucos^, where C is given by 

cos C = sin a sin /3 , (9) 

in terms of a and the angle /3 between k and the direction to the star's spin pole projected onto the 
plane tangent to the stellar surface at the point of the emission. 
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In order to compute the waveform seen by a distant observer located at an inclination ^obs relative 
to the star's spin axis, we divide the hot spot on the stellar surface into a fine grid in colatitude and 
longitude. We then consider the flux from an infinitesimal emitting area dA'^ in the comoving frame 
around each grid point. We transform this flux into the flux in the static frame using equation (j5]) and 
the relation 

dAi = SdA'i (10) 

between the infinitesimal area dA'- measured in the comoving frame at the stellar surface and the in- 
finitesimal area dAi measured in the static frame at the stellar surface. Next we use spherical trigonom- 
etry and ray-tracing in the Schwarzschild spacetime to determine the direction of emission (aj,/3j), 
measured in the static frame at the stellar surface, that is required for a light ray originating at a given 
grid point to reach the observer. We determine the angular separation ipi between the location {9i,ipi) 
of the grid point and the direction to the observer using the spherical trigonometric relation 



cos ipi = cos 9i cos ^obs + sin 9i sin 6'obs cos ifi 



and then use the implicit relation 



sin a,- dx 



y(r^-2W^r^T^^^2M^7Ry^2^h?^ 



(11) 



(12) 



between ipi and Oi to determine aj. We accurately evaluate the integral in equation (fT2]) using a 
combination of analytical and numerical methods (see Section lA.l.ip . Finally, we determine /3j using 
the spherical trigonometric relation 



cos f3i = (cos ^obs — cos 9i cos Ipi) /{sin 9i sin ipi) . 



(13) 



Given the position of the observer and the emitting point on the star, this algorithm allows us to solve 
in a single step for the angles Oj and /3i at which the emitted ray leaves the stellar surface. 

Knowing ai and /3j, we can evaluate the energy-resolved photo n number flux arriving at the distant 
observer from the inflnitesimal element dAi, using the expression (jPoutanen fc Gierlinskill2003l . eq. 13) 



dFiiE) 



[l + z) 



^iIo{Eo,ai) dAi cos ai 



E 



1)2 



sin ^ dip 



sin a da 



-1 



(14) 



where E = {1 + z) ^£"0 is the photon energy measured in the static frame at infinity and 1 -\- z = 
(1 — 2M/i?)~^'2 is iiiQ gravitational redshift from the stellar surface to infinity. 

In order to compute efficiently the waveform produced by the emission from the entire hot spot, 
we proceed in two steps. First, we determine each of the waveforms produced by emission from a set 
of emitting areas {d^'j} with tp = that span the colatitudes within the hot spot. In doing this, we 
take into account the fact that light emitted from different locations {R, 9i,ipi) takes different lengths of 
time to reach the observer. However, the waveform seen by the observer depends only on the difference 
in the light travel time from different emitting points. (A different choice for the zero of time or a 
different total propagation time would shift the arrival time of waveform but would not affect the shape 
of the waveform, which is what concerns us here.) Hence, we need to compute only the differences 
in the arrival times of the different rays that reach the observer. We choose to compute the arrival 
time of each ray relative to the arrival time of a radial ray emitted from a point on the stellar surface 
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immediately below the observer, which in our time coordinate is t = (see above). The arrival time of 
a photon that leaves the surface at an angle a. to the normal as measured in the static frame is then 



Af= dr - 

Jo x2 (1 - 2Mx/R) 



1 - sin^ a (1 - 2M/R)~^ (1 - 2Mx/R) x"^ 



(15) 



We accurately evaluate this integral using a combination of analytical and numerical methods, as de- 
scribed in Appendix El Once we have determined the waveform produced by emitting areas at (/? = 
that span the colatitudes within the hot spot, we compute the full waveform by adding the waveforms 
produced by emission from the grid points that have different values of tp. These waveforms can be 
generated quickly by appropriately shifting the phase of the waveform produced by the corresponding 
emitting element at 99 = 0. 

We have validated our waveform code using a suite of code tests. These tests, and the results, are 
discussed in Appendix |Al We have carefully chosen values for the integration, hot spot, and angular 
resolution parameters that we use in our waveform code to provide resolutions fine enough to meet our 
accuracy requirements, but no finer, so that our code runs as fast as possible. 



3.2. Bayesian Analysis and Sampling Methods 

We wish to determine both the most probable ("best-fit") values of the parameters in our waveform 
model, given an observed burst oscillation waveform, and the confidence regions for the values of these 
parameters. Both goals can be accomplished simultaneously and efficiently using Bayesian inference 
and an MCMC algorithm to sample the parameter space. Once we have computed the best-fit values 
of the model parameters, we can determine the accuracy of the fit by comparing them with the values 
that were used to generate the synthetic observed waveform. 

The most probable values of the parameters y in our waveform model and their confidence regions 
can be determined using the posterior probability distribution p{y\D,I), where D is the synthetic 
energy- and oscillation phase-resolved waveform data of interest and I is any information obtained 
prior to the measurements under consideration. The desired posterior probability distribution can be 
obtained from the likelihood of the data, given the parameter values, using Bayes' theorem 

p{y\DJ)<^p{D\y,I)p{y\I), (16) 

where p{y\I) is the prior probability distribution of the parameter values and the constant of propor- 
tionality is the inverse of the normalization factor. This constant of proportionality is irrelevant when 
estimating the values of the parameters in a given model. In the present analysis, we use the most 
uninformative prior, namely, we assume that p(y|/) is uniform for parameter values within the physical 
ranges we consider. Then, for parameter values within these ranges, 

p{y\DJ)^p{D\yJ). (17) 

In performing the MCMC sampling of the parameter space, we use first-order Markov chains and 
start from a random point in the parameter space. At each step, we generate a proposed new set of 
parameter values y' based on the current set of parameter values y'") by drawing from a proposal 
distribution Pp(y'|y^"^). The proposed new set of parameter values is then accepted or rejected with 
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specified probabilities. We use the basic Metropolis algorithm (see, e.g.. Ivon Toussaint]|201ll ). drawing 
the new set of parameter values from a joint-normal distribution 

Pp(y |y^"^)=Pp(y^"^|y')~A^(y -y^"^'^) , (is) 

where the elements of a are the standard deviations for each parameter, to be specified. We adopt the 
acceptance probability (which is also the transition probability) 

With this choice, we always accept the proposed new set of parameter values if its probability is 
higher than that of the current set; otherwise, we accept the proposed new set with probability 
p(y'|-D,I)/p(y("^|-D,/). As required, this transition probability satisfies the detailed balance condition, 
as can be readily verified. Using relation (J17p . this ratio can be written 

v{y'\DJ) piD\y',I) ^2^^ 



p{y(^)\D,l) p{D\y(^),l) ' 

Hence, at each step the acceptance probability for the proposed set of parameter values is determined by 
calculating the ratio of the likelihoods of the data given the current parameters and given the proposed 
parameter values. 

If Poisson noise is the only source of fluctuations in the data, the likelihood of the "observed" data, 
given a particular set y of values for the model parameters, is 

C^p{D\y,I) = ll^^^^e--^(y\ (21) 

i 

where the product is over all the oscillation phase and energy bins, di is the measured number of counts 
in the i^^ bin, and rn-j(y) is the number of counts in the i^^ bin predicted by the model for the trial set 
y of parameter values. 

In our MCMC algorithm, we use the log likelihood only to determine the transition probability, 
which depends only on the difference of log likelihoods, so the dj! terms in the log likelihood cancel out, 
producing the expression 

log £ = ^ di log mi{y) - A^modcily) , (22) 

i 

where 

iVmodcl(y)= J^rni(y) (23) 

i 

is the total number of photon counts in the model spectrum. Unlike many situations, where the 
normalization of the model is a free parameter and hence the total number of counts can be adjusted 
to be the same for every set of parameter values, here the total number of counts, and hence the 
normalization, depends explicitly on the distance to the star, the angular diameter of the spot, and 
other model parameters. The normalization of the model is therefore a key quantity for discriminating 
between different sets of parameter values. 
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3.3. Estimating M and R 

We wish to determine the best-fit values and confidence regions of the parameters M and R in our 
waveform model, given a synthetic observed waveform of interest. With a uniform distribution over 
the allowed values of the model parameters as our prior, the posterior probability of a particular set 
of parameter values y, given the data D, is proportional to p{D\y), the likelihood of the data given 
those values of the parameters. There are two main computational tasks: (1) computing the relative 
likelihoods of the data over a set of trial model waveforms chosen to span and adequately sample the 
parameter space of the model, and (2) marginalizing the resulting posterior distribution by integrating 
over all the parameters in the model except M and R. 



3.3.1. Construction of synthetic observed waveforms 

The data D representing a synthetic observed waveform is a list di of the (integer) number of counts 
observed in each oscillation phase-energy bin. We produce each synthetic observed waveform in three 
steps. First, we use the waveform generating code described in Section [3T] to compute the oscillation- 
phase- and energy-resolved waveform for a set of model parameter values of interest. Second, we add 
phase-independent (but energy-dependent) counts from our background model, which as we explained 
in Section 12.11 is a catch-all intended to mimic possible emission from the entire stellar surface, an 
accretion disk, and other sources in the field of view, as well as the instrumental background. Finally, 
we Poisson-sample the total number of photons in each of the phase-energy bins. 



3.3.2. The computational problem 

The straightforward way to determine the joint posterior probability distribution of M and R would 
be to integrate the joint distribution of all the parameters in our waveform model over every parameter 
except M and R. This integral. 



p{M, 



R\D,I) = j AV p{y\DJ), (24) 



could in principle be computed using a Monte Carlo algorithm. Here V is the volume of the model 
parameter space when the M-R subspace is excluded. If we were to sample the integrand at A'^ points 
{xi} picked randomly but uniformly within V, the uncertainty in p{M,R\D.,I) would be (Press et al. 
1999, Section 7.6) 

Ap(M, R\D, I) « y^MZM! , (25) 

where 

1 ^ 

j=0 
and 

1 ^ 
{p')^j^Y.P^x,\Djf. (27) 

i=0 
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Expression (|25p shows that the uncertainty in the marginahzed posterior probabihty distribution ob- 
tained by Monte Carlo integration over the fuh posterior probability distribution decreases as the number 
A^ of sample points increases, but only as 1/\/N . 



3.3.3. Computational procedure 

In this work we seek to determine the most probable values of M and R and their confidence 
intervals by comparing our waveform model with a synthetic observed waveform, using a Bayesian 
approach. Each model waveform is a list of the expected number of photons rrii in each oscillation 
phase-energy bin. We construct a complete model waveform by computing a hot spot waveform (i.e., 
a model waveform without any background) and then adding a model of the background. For the 
waveform model we use here, specifying a complete waveform requires specifying 38 model parameters: 
M, R, the triplet of angles ^spoti ^^spoti and ^obs (which define the parameter subspace y'), plus 
the color temperature T^o of the emission at the stellar surface measured in the comoving frame, the 
distance d to the star, the absolute phase of the oscillation (J)q, and the background counts in 30 energy 
channels. Determining the joint posterior probability distribution of M and R using equation (I24p 
therefore requires accurate computation of the posterior probability distribution over a high-dimensional 
parameter space and subsequent computation of the marginalization integral over this space. We found 
that the computational effort needed to achieve sufficient accuracy using this approach was excessive. 
We therefore sought a more efficient approach, which we now describe. We presume that when the 
data from a large-area timing mission become available, the computational resources needed to do full 
Bayesian analyses of these data will also be available. 

The purpose of our analysis here is not to reproduce all the steps that would be needed for a full 
Bayesian analysis of a real observed waveform, but rather to determine the precision and accuracy with 
which such a full analysis could determine M and R. Our initial exploration of the computational 
problem revealed several shortcuts that we could use for the current study that substantially reduce 
the computational burden but do not alter the results significantly. Using these shortcuts, the most 
probable values of M and R and their confidence intervals can usually be computed for a single synthetic 
observed waveform in 50-100 clock hours, running a parallel code on 150 nodes of a fast CPU cluster. 

The procedure we use to determine the most probable values of M and R and their confidence 
intervals has seven steps: (1) construct an initial grid of points in the M-R plane; (2) for each M-R pair 
in this grid, choose values for the spot inclination ^spot) the spot angular radius A^spoti and the observer 
inclination ^obsj (3) determine the color temperature T^o of the emission from the hot spot measured in 
the comoving frame at the stellar surface, the absolute phase (j)o of the waveform, the distance d to the 
star, and the background model that maximize the likelihood of the observed waveform; (4) sample the 
likelihood distribution p[D\y') over the three-dimensional parameter space y' consisting of the three 
angles ^spot, A^spot, and ^obs! using the most probable values of Tco, (/'o, d, and the background for each 
angle triplet; (5) integrate the posterior probability distribution p(y|Z?) over y', using the most probable 
values of the other model parameters except M and R, and thus associate an integrated probability 
density p{M,R\D) with each point in the M-R grid; (6) refine and extend the M-R grid as needed 
to determine the most probable values of M and R and their confidence intervals with the desired 
accuracy; (7) use this final integrated probability density p{M,R\D) to determine the most probable 
values of M and R and their la and 3cr Bayesian confidence intervals. We now explain several aspects 
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of this procedure in more detail. 

Sampling M and R. To compute the posterior probabihty distribution p{y\D) efficiently, we first 
construct an initial grid of points in the M-R subspace surrounding the values of M and R that were 
used to generate the synthetic observed waveform being analyzed. This approach will of course not be 
possible when analyzing real data, because the actual values of M and R will not be known in advance. 
Instead, one will have to search the entire M-R parameter space, greatly increasing the computational 
effort required to adequately sample p{y\D). Once we have computed p{y\D) on the initial M-R grid, 
we use the results to extend and refine the grid, repeatedly if necessary, until the probability density 
has been determined over the relevant portion of the M-R plane with sufficient accuracy that we can 
determine the best-fit values of M and R and their uncertainties with the desired accuracy. 

Determining the emission color temperature T^o- Our preliminary analysis showed that T^o is 
strongly constrained by the spectrum of the phase-dependent part of the observed waveform, and is 
much more strongly constrained than the angles ^spot, A^spot, and ^obs- Hence the range of Tco values 
where the probability density is appreciable is very small and using the MCMC algorithm to sample 
this parameter is therefore very inefficient. Consequently we used a different approach to determine the 
value of Tco in our trial waveforms. 

All the trial waveforms used in this paper were generated assuming that the spectrum of the 
emission from the hot spot has the same shape as the Planck spectrum. Most of the synthetic observed 
waveforms that we consider here were generated assuming that the spectrum of the emission from the 
hot spot has the same shape as a Planck spectrum. However, a few of the synthetic observed waveforms 
considered here were generated assuming a hot spot emission spectrum having the same shape as a 
Bose-Einstein spectrum, to explore the effects on the estimates of M and R of fitting a model that 
assumes a spectrum somewhat different from actual spectrum. We use different approaches to generate 
trial values of Tco in these two cases. 

When the comoving emission spectrum used to generate the synthetic waveform has a Planck shape, 
we found that using the redshift relation for a non-rotating star, namely, 

T,, = T^{l-2M/R)-^/\ (28) 

to estimate Tco using the observed radiation temperature Too gives Tco to within 5%, for spot rotation 
frequencies < 600 Hz. Consequently, to reduce the computational burden during our large-scale runs, 
for this synthetic spectrum we assume relation (I28p and use as our trial value of Tco the value given by 

_ (1 - 2Mo/i?o)V2 

^co-ico,0 (i_2M/i?)l/2 ' ^^^^ 

where the quantities with '0' subscripts have the values that were used in generating the synthetic 
waveform. Using this approach decreases the time required to compute p{M,R\D) by a factor of five 
compared to the time required to solve for the most probable value of Tco. 

When the emission spectrum used to generate the synthetic waveform has a Bose-Einstein shape, 
we included Tco as one of the parameters determined by maximizing the likelihood of the observed 
waveform, as we discuss below. 

Changing the blackbody spectrum in the model waveform takes very little computational effort. 
The reason is that changing from one Planck spectrum with a given temperature and normalization to 
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another with a different temperature and normahzation can be done simply by mapping the spectrum 
to the different energy and multiplying it by the factor needed to renormalize it. Because photons 
of all energies are Doppler shifted, gravitationally redshifted, and deflected in the same way, once the 
energy-dependent waveform for one comoving emission temperature has been computed, the waveform 
for any other comoving emission temperature can be computed using the appropriate mapping and 
renormalization factor. Thus, the very time-consuming step of tracing rays needs to be done only once. 
To preserve accuracy when performing the mapping and renormalization, it is important to use tables 
that span the required energy range, so that interpolation can be used and extrapolation is not required. 

If there is an atomic scattering line in the synthetic hot spot spectrum, we assume that we know 
the rest energy, width, and optical depth of the line and use this information to include a line with 
these properties in the spectrum of the waveform model. 

Determining the most probable values ofTco, d, (J)q, and the background. Given an M-R pair, there 
are 35 additional waveform parameters, if the synthetic emission has the shape of a Planck spectrum 
and Tco is determined as described above, or 36 additional parameters, if the synthetic emission has the 
shape of a Bose-Einstein spectrum and T^o has to be determined by a likelihood analysis. 

As noted above, our preliminary analysis showed that Tco is strongly constrained by the spectrum 
of the phase-dependent part of the observed waveform. Hence, if Tco must be determined, we do so by 
maximizing the likelihood of the observed waveform, given the model waveform, using bisection. 

In our preliminary investigation, we found that when the other waveform parameters are held fixed, 
the distance d to the star (which is equivalent to the normalization factor of the Planck spectral shape) 
is also very tightly constrained by the total number of observed counts. The total number of observed 
counts is in turn very tightly constrained by the data. Consequently, determining d by MCMC sampling 
and marginalization is very inefficient. We therefore determine d by maximizing the likelihood of the 
data using bisection. This method determines the most probable value of the distance to within 1%. 

Our preliminary analysis showed that the phase (pQ of the oscillation is very tightly constrained by 
the observed waveform. We therefore determine (J)q by maximizing the likelihood using bisection. The 
first step is to apply the same, tentative phase shift to the model waveform in each photon-energy bin. 
We do this by Fourier transforming the trial waveform at every photon energy, applying the tentative 
phase shift to all the harmonic components, and then inverting the Fourier transform to obtain the 
shifted model waveform. We perform the phase shift in the frequency domain rather than in the time 
domain to increase the accuracy: a phase shift applied in the time domain would have to be at least 
as large as the frequency times the width of a time bin, whereas even a very small phase shift can be 
applied to every harmonic in the frequency domain. Then, when the Fourier-transform is inverted the 
full waveform will be shifted by this small amount. 

We compute the background model using bisection to determine the number of background counts 
at each energy that, when added to the phase-shifted hot spot waveform, produces a complete model 
waveform that maximizes the likelihood of the observed waveform. 

Sampling the three angles. The values of the three angles 6'spot5 A^gpot, and ^obs are much less 
tightly constrained than are Tco, d^ (po, and the background. Therefore, in our initial sampling of these 
angles we use the MCMC algorithm outlined above. In this step, we typically generate 30 chains, each 
with a length of 40, for a total of 1200 triplets. For each triplet of angles, we simultaneously determine 
the values of the remaining 32 (or 33) parameters that maximize the likelihood, using bisection. We 
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then select the triplets that have a log likelihood within 20 of the maximum log likelihood. Finally, we 
use this collection of triplets to define the ranges of the three angles that we sample in the next step of 
our procedure. 

For our second, final sample of the three angles, we typically choose 10,000 angle triplets randomly 
and uniformly within the ranges determined in the previous step. The number of triplets was chosen to 
provide enough resolution to obtain an acceptably accurate value of the integrated probability density. 
For each triplet of angles, we again determine the values of the remaining 32 (or 33) parameters that 
maximize the likelihood using bisection. We then sum the probabilities at all these triplets, to associate 
an integrated probability density p{M,R\D) with this particular M-R pair. 

Completing the sampling of M and R. Having computed the integrated probability density for 
a particular M-R pair, we then choose a different M-R pair in our grid and repeat the process. If 
necessary, we extend and refine the M-R grid, in order to determine the most probable values of M and 
R and their la and 3(T confidence intervals with the desired accuracy. We use a common normalization 
for the integrated probability density of different M-R pairs, so it can be compared for different pairs. 
Finally, we use the resulting probability distribution in the M-R subspace to determine the most 
probable values of M and R and their la and 3(T confidence intervals. We tested this algorithm in 
numerous cases and found that it does a good job of determining the most probable M-R pair and the 
corresponding confidence regions. 
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4. WAVEFORM ANALYSIS AND RESULTS 

Our burst oscillation waveform analysis proceeds in three steps: 

1. Generate a synthetic observed waveform, using the method described in Section [3.11 

2. Compute the joint posterior distribution of the mass and radius parameters in our waveform 
model, given this synthetic waveform, using the computational methods described in Sections 
and 



3. Use this posterior distribution to determine the most probable values of M and R and the la and 
3(T confidence regions in the M-R plane, given this synthetic waveform. 

We use this procedure to explore how the accuracy and precision of M and R estimates depend on the 
physical characteristics of the observed system by generating a number of different synthetic waveforms 
corresponding to systems with different physical characteristics, such as spot rotation rate, spot incli- 
nation, and viewer inclination. We then explore the effects on the uncertainties in M and R of having 
additional information about the observed system that is independent of that provided by the burst 
oscillation waveform. We consider knowledge of the distance, the observer's inclination, and the size and 
spectrum of the background, measurement of the properties of an atomic scattering line in the emission 
from the hot spot, and high-precision measurements and modeling of the spectrum of the X-ray bursts 
from the star in question. We also explore the effects on M and R estimates if the actual spot shapes. 
X-ray spectra, and radiation beaming patterns differ from those assumed in our waveform model. Fi- 
nally, we investigate the constraints that can be obtained by analyzing data on burst oscillations with 
changing waveforms, collected during a single burst or during multiple bursts. 



4.1. Waveform Analysis 

All the synthetic waveforms used in this work assume M = 1.6 Mq and R = 5.0 M. We compute 
the joint posterior distribution of the parameters M and R in our waveform model, given the synthetic 
waveform, for 3.5 < R/M < 6.8. Given that all the synthetic waveforms assumed R/M = 5.0, this 
range in R/M is adequate to evaluate all cases that provide relatively tight constraints on M and R. 



4.1.1. Model burst oscillation waveforms 

The hot spot emission model that we usually use in fitting the synthetic waveform data assumes 
that the spot is a uniformly emitting circular area on the stellar surface that is rotating around the 
star and that each point on the surface of the hot spot radiates with a spectrum that has the shape of 
a Planck spectrum and the angular pattern expected for emission from a semi-infinite, nonrelativistic, 
electron scattering atmosphere (see Section [3T] for further details). 

The full waveform model has 9 parameters that describe the waveform produced by the emission 
from the hot spot (the star's gravitational mass M, circumferential radius R, and the spot rotation 
frequency Umt', the spot inclination or colatitude 6'spot and angular radius A6'spot; the color temperature 
Tco of the spot emission as measured in the comoving frame at the stellar surface; the observer's 
inclination ^o^s relative to the rotation axis and the distance d to the star; and the phase of the 
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waveform). The waveform model has an additional 30 parameters that describe the background model 
(the number of background counts in each of 30 energy channels). Thus there are in principle a total 
of 39 parameters in the model. 

In practice, the spot rotation frequency can be determined separately and much more accurately 
than any of the other parameters in the waveform model. Therefore, to speed up the computations 
for the present study we assume that i/fot is already known. As discussed in Section 13.3.31 ^co can be 
determined very accurately using the color temperature observed at infinity when the spectrum of the 
emission from the hot spot that produces the observed waveform has the shape of a Planck spectrum. 
As discussed there, in this case we use this fact to determine the value of Tco separately from the fitting 
procedure, to speed up the computations by a factor of five. When the spectrum from the hot spot 
does not have the shape of a Planck spectrum, Tco must be included as one of the fitted parameters. 
Thus, the hot spot component of the waveform model has 7 adjustable parameters, if the spectrum has 
the shape of a Planck spectrum, or 8 free parameters, if it does not. The full waveform model, which 
includes the background, therefore has either 37 or 38 adjustable parameters. Each synthetic oscillation 
waveform has 480 data points (see below), so the number of degrees of freedom in the fits described 
here is either 480 — 37 = 443 or 480 — 38 = 442, depending on the hot spot spectrum that was assumed 
in generating the synthetic waveform data. 

If the properties of the observed neutron star are consistent with the properties assumed here, we 
expect the constraints on M and R to improve as the ratio TZ of the number of modulated counts to 
the fluctuations in the total number of counts (see equation ([T])) increases. For a fixed value of 7^, 
we expect the constraints to be tighter when the spot and observer are closer to the star's rotational 
equator, because the relativistic Doppler shift and aberration are then greater, producing a waveform 
that depends more sensitively on the radius of the star. 



4-1-2. Synthetic burst oscillation waveforms 

Each synthetic observed waveform is represented by the number di of counts in each of 30 energy 
bins of width 0.3 keV covering 3.5 keV to 12.5 keV at each of 16 rotational phases. There are thus a 
total of 480 data points in each synthetic waveform. 

The physical models we use to generate the component of these waveforms that is produced by the 
hot spot are characterized by the gravitational mass M and circumferential radius R of the star; the spot 
rotation frequency frot! the location, size, and shape of the hot spot; the spectrum and beaming pattern 
of the emission from the hot spot; and the distance d to the observer and the observer's inclination 
^obs relative to the star's rotation axis. The waveforms used in this work assume that the emission 
from the hot spot has the shape of a Planck spectrum or a Bose-Einstein spectrum with a temperature 
kTco = 2.0 keV, as measured in a frame comoving with the surface of the rotating hot spot. To facilitate 
comparison of different cases, we adjust the distance to the star so that the expected number of counts 
from the hot spot is 10^ in all cases. 

To explore the effects of background counts, we add a background component to the waveform 
produced by the hot spot. As discussed in Section [2TT1 this component is a catch-all for all counts not 
produced by radiation from the hot spot. These counts could be produced by emission from unassociated 
sources in the field, the accretion disk, the non-spot portion of the star, instrumental backgrounds, or 
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any combination of these. We assume that the background does not vary at frequencies commensurate 
with the spot rotation frequency. We model the background by assuming uniform emission from the 
entire surface of the star with a spectrum having the shape of a Planck spectrum with a temperature of 
1.5 keV as seen in a frame comoving with the surface of the star, which for this purpose is assumed to 
be rotating with the same frequency as the hot spot. To study the effects of different background levels, 
we adjust the normalization of the background component so that the number of counts expected from 
the background is either 0.3, 1.0, or 9.0 times the 10^ counts expected from the hot spot. We refer to 
these as our low, medium, and high backgrounds. 

Once the preliminary synthetic waveform is generated, it is Poisson sampled to produce the complete 
synthetic waveform. Because the waveforms are Poisson-sampled, the actual counts contributed to the 
complete waveform by the hot spot and the background are slightly different from their expected values. 

Standard synthetic waveforms. We first analyze synthetic waveforms generated using burst oscilla- 
tion models with the same qualitative properties as the model that we fit to the data. These "standard" 
synthetic waveforms assume that the hot spot is circular and uniformly emitting, that the spectrum of 
the emission from the hot spot has the shape of a Planck spectrum, as measured in the comoving frame 
at the stellar surface, and that the beaming pattern of the emission from the hot spot is that expected 
for a nonrelativistic election scattering atmosphere, which is described by the Hopf beaming function 
(see Section 13. ip . The standard waveforms we analyze here assume that the angular radius A^gpot of 
the hot spot is 25°. 

Reference cases. The hot spot and observer inclinations have the biggest effects on the uncertainties 
in estimates of M and R. To facilitate comparisons, we constructed several standard waveforms for a 
"low-inclination" reference case and for a "high-inclination" reference case. The model parameters that 
define these two cases are frotj ^obs; and ^spot- The values of these parameters for the two reference 
cases are listed in Table [H 

The low-inclination reference case is motivated by the evidence that the m agnetic poles of ma ny 



accreting neutron stars with millisecond spin periods are close their spin poles (see lLamb et al.ll2009s 
This evidence includes the low fractional amplitudes of the X-ray oscillations of accretion-powered 
millisecond X-ray pulsars, their highly sinusoidal waveforms, and their intermittency in some stars. In 
the low-inclination reference case, the hot spot and the observer are at inclinations of 20° and 60°, 
respectively, and the rotation frequency is 400 Hz. 

Table 1. System parameters for the two reference cases'^ 



Case label i^rot (Hz) 6'obs (deg) 6'spot (deg) 

low inclination 400 60 20 

high inclination 600 90 90 

'^Here t'rot is the rotation frequency of the hot spot, ^obs 
is the inclination of the observer relative to the rotation 
axis, and ^spot is the inclination (colatitude) of the center 
of the hot spot. 
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Even if the star's magnetic poles are close to its spin pole, the region heated by nuclear burning or 
surface modes may not be at or even near a magnetic pole, and hence may not be close to the spin pole. 
As an example, the very high-amplitude X-ray oscillations found near the beginning of a burst from 
4U 1636—536 bv lStrohmaver et al.1 (jl998l ) are difficult to understand if the heated region is near the spin 
pole. We therefore also consider a high-inclination reference case in which the hot spot and observer 
inclinations are both 90° (i.e., both are in the rotational equator) and the spot rotation frequency is 
600 Hz. These choices produce the largest expected Doppler shifts and relativistic aberration, and 
should therefore illustrate the tightest possible constraints on M and R. 

Nonstandard synthetic waveforms. In addition to our "standard" synthetic waveforms, which as- 
sume a circular hot spot with a 25° angular radius, a Planck emission spectrum, and the radiation 
beaming pattern appropriate for emission from a Thomson scattering atmosphere (described by the 
Hopf function), we also generated several other waveforms, to explore the consequences of fitting an 
incorrect waveform model to the data. These "nonstandard" waveforms were generated for spots elon- 
gated in the north-south and east-west directions, for a Bose-Einstein emission spectrum, or for isotropic 
beaming from the surface of the hot spot. 



4-. 1.3. Independent additional information 

We explore the improvements in constraints on M and R that become possible with independent 
additional information about important system properties. There are good prospects for determining 
the angle of our line o f sight relative to the a xis of the binary o r bit using optical burst reflection 
mapping ( Casaresll20ld ) and Fe-line modeling ( Cackett et al.ll2010l : lEgron et al.ll201ll ). The neutron 
stars in bursting systems ar e thought to have acquired nearly all their angular momentum via accretion 
(JLamb fc BoutloukosI l2008l ) . so the spin axes of these stars are likely to be closely aligned with the 
orbital axes of their systems, allowing us to infer the inclination of the observer relative to the stellar 
spin axis. The distance of a bursting neutron star is often tricky to measure, but it can be constrained if 
the system is in a g lobular cluster. P arallax distances for such stars are likely to become available in the 
future from GAIA (JEver et al.ll2013l ). The properties of the non-spot background are likely to be difficult 
to constrain independently of the waveform fitting process, but if the spectrum of the background is 
sufficiently distinct from the spectrum of the spot radiation (e.g., if both have Planck spectral shapes but 
have easily distinguishable temperatures), it may be possible to restrict the properties of the background 
using spectral measurements. It is also possible that an atomic resonance scattering line could be 
identified in the emission from the hot spot. We expect that having information like this would improve 
the constraints on M and R. 

It may be possible to obtain other independent information by fitting high-precision spectral models 
to the observed spectra of long bursts or jointly fitting such models to the spectra of many bursts from 
the same source. At high burst temperatures, the burst atmosphere is expected to be scattering- 
dominated and to produce spectra with shapes that are very close to the shape of a Bose-Einstein 
spectru m (Lo et al., in prepar ation). This expectation is consistent with the observed spectra of such 
bursts (JBoutloukos et al.ll2010l ). If the observed spectra had exactly a Bose-Einstein spectral shape, they 
would provide no information about the properties of the star, because a thermodynamic equilibrium 
spectrum provides no such information. However, we have found that detailed, high-precision model 
atmosphere spectra provide a much better description of the best available spectral data (the RXTE 
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PCA data on the superburst from 4U 1820—30) than do Bose-Einstein spectra, so some information 
about the properties of the star can be obtained. Information about the mass and radius of the star 
comes mainly from the smah deviation of the spectrum from a Bose-Einstein shape at low photon 
energies that is caused by the free- free opacity. This deviation allows one to determine (1 -|- z)/g^'^ , 
where z is the surface redshift and g is the surface gravity of the star. This relation between M and R 
complements the constraints on M and R obtained from fits to the burst oscillation waveform. 



4- 1.4- Errors in the waveform model 

Most of the results we present here assume that we know the shape of the hot spot and the spectrum 
and beaming pattern of the radiation from the spot. However, in practice these properties may not be 
accurately known. If any important properties of the actual star differ significantly from the properties 
assumed in the model being fit, we are attempting to fit the data using an incorrect model. In this 
situation we expect the fit to be poor and the best-fit values of M and R to differ from their actual 
values. Of particular concern is the possibility of situations in which the fit is formally good but the 
best-fit values of M and R are significantly different from their actual values. A full analysis of possible 
systematic errors is beyond the scope of a this work, but we briefiy explore the effects of errors in the 
model of the spot emission. 



4.2. Constraints on M and R 

Consider now the constraints that can be obtained by fitting a waveform model to an observed 
burst oscillation waveform. Our results show that the uncertainties in M and R estimates produced by 
statistical fiuctuations in the waveform are typically only a percent or two, if ~ 10^ counts are collected 
from the hot spot of a given bursting neutron star and the values of all other parameters are known 
independently of the waveform analysis. However, for realistic situations in which some of the other 
parameters in the model must be determined from the waveform, the statistical uncertainties are usually 
much larger. The reason for this is that the effects on the waveform of changing different parameters 
in the waveform model are often very similar. These degeneracies of the waveform with respect to 
changes in the values of the model parameters are an inherent property of any physical model based on 
a rotating hot spot and cannot be removed by "improving" the model. These degeneracies infiate the 
uncertainties in M and R estimates produced by statistical fiuctuations. 

The synthetic waveforms we consider here and the resulting la uncertainties in M and R are 
summarized in Table [2j Plots of the la and 3a uncertainty contours for each of these waveforms are 
presented and discussed below. The x^/dof for the best-fit model is shown on each plot. 

The uncertainties in M and R quoted in this section were estimated by projecting the extent of the 
two-dimensional uncertainty regions onto the M and R axes. This approach neglects the fact that the 
probability distributions of M and R are often correlated. It also somewhat overestimates the individual 
uncertainties in M and R, because the central peak in the distribution of the unwanted quantity would 
centrally weight the probability distribution of the quantity being considered, if the unwanted quantity 
were actually marginalized. This effect is neglected when the uncertainty region is simply projected. 
We note that the quoted uncertainties are themselves slightly uncertain, due to incomplete sampling of 
the high-dimensional parameter space. 
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Table 2. Synthetic waveforms, fitted models, and corresponding figures'^ 



Figure Inclinations Background'^ Additional description 
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^spot = 80°, all parameters unknown 

^spot = 60°, all parameters unknown 

observer inclination known 
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atomic line present in the data 

atomic line present in the data 

background known 

background known 

hot spot elongated longitudinally 

hot spot elongated latitudinally 

Bose-Einstein spectrum in the data 
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^AU synthetic waveforms were normalized to produce ~ 10^ counts from the hot spot. ^See Table [T] for 
the parameters used in generating the high- and low-inclination synthetic waveforms. The elongated spots 
are defined in the text. '^The backgrounds in all synthetic waveforms are phase- independent, with the 
spectral shape of a Planck spectrum having a temperature at infinity of kT = 1.5 keV. Our low, medium, 
and high backgrounds produce ~ 3 x 10^, ~ 1 x 10^, and ~ 9 x 10^ detected counts. Approximate la 
fractional uncertainties in M and R estimated by projecting the la contours onto the M and R axes (see 
text). '^These purely statistical uncertainties assume that all relevant properties of the actual system are 
accurately described by the waveform model, that the values of these properties lie within the range of 
values covered by models considered in the waveform fitting process, and that the values of all the model 
parameters except M and R are known exactly by some means that is independent of the waveform 
fitting process, thereby eliminating the degeneracies of the model parameters. -^Reliable estimates of the 
uncertainties are not possible because the confidence regions are too large. ^The reduced x^ of this fit 
indicates that the model is incorrect. 
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Fig. 1. — These contours are based on the highly unreahstic assumptions that ah the relevant properties 
of the actual system are correctly included in the waveform model, that the actual values of these 
properties lie within the range of values possible in the waveform model, and that the values of all the 
parameters in the waveform model except M and R are known exactly, independent of the waveform 
fitting process. This procedure eliminates completely the effects of the degeneracy of the waveform 
with respect to variations in the values of the system parameters. The + symbol indicates the radius 
and mass that were used in generating the synthetic waveform and the star symbol shows where the 
marginalized posterior probability density is highest. The dotted and solid curves show, respectively, 
the 1(7 and 3a confidence contours. The reduced x^ indicates that the fit is acceptable. See the text for 
further details. 



4-2.1. Effects of statistical fluctuations 

We can estimate the uncertainties in estimates of M and R produced by the statistical fluctuations 
in the observed waveform by considering the (highly unrealistic) situation in which all the relevant 
properties of the actual system are correctly included in the waveform model, the actual values of these 
properties lie within the range of values possible in the model, and the values of all the parameters in 
the model except M and R are known exactly by some means that is completely independent of the 
waveform fitting process. 

Figure [1] shows an example in which the assumptions in the model used to fit the synthetic waveform 
data are identical to the assumptions in the model that was used to produce the waveform data. The 
waveform data were generated using the parameter values of the high-inclination reference case and 
our medium background. In fitting the model to the waveform data, the values of all the parameters 
in the model except M and R were set to the values used to produce the synthetic waveform. This 
procedure completely eliminates the effects of the degeneracy of the waveform under variations of the 
model parameters. The most probable values values of M and R were then determined using our 
standard fitting and marginalization procedures. The la uncertainties in M and R are about 1.9% and 
1.3%, respectively, whereas the 3a uncertainties in M and R are about 5% and 3%. The value of x^/dof 
indicates that the fit is acceptable. 



These results can be compared with those of IStrohmayed (|2004l ). who investigated the constraints 
on M and R that could be obtained by fitting a bolometric waveform model to synthetic bolometric 
waveform data. He assumed that the hot spot and the observer are both in the rotation equator 
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(their most favorable locations) and that the size, shape, and location of the hot spot, the spectrum 
and beaming pattern of the radiation from the hot spot, and the inclination of the observer are all 
known independently of t he waveform fitting analysis, so that only M and R have to be estimated 



from the waveform data. IStrohmaveii a lso assumed there is no background. With these assumptions 
and approximations, IStrohmaveiJ (|2004l ) found that fitting a bolometric waveform model to bolometric 
waveform data from an X-ray instrument with an effective area ~ 10 m^ could provide estimates of M 
and R with uncertainties of a few percent. 

Our simulations show that the statistical uncertainties in M and R scale roughly as 1/7^, where 
TZ = A'^Qsc/v'^Vtot in terms of the number Nqsc of counts in the oscillating component of the synthetic 
waveform and the total number A'tot of counts in the waveform (see equation ([1]) and the discussion 
following it). This is illustrated in Table [3l In these particular examples, the expected number A^back of 
counts in the background is equal to the number A'^gpot of counts from the hot spot, so A^osc oc A'spot oc 
A^'spot + A^'back = A^'tot and hence TZ is proportional to i/A'tot- 



4-2.2. Effects of parameter degeneracies 

Comparison of the uncertainties listed in Table [2] and comparison of the confidence regions shown 
in Figure [1] with the confidence regions shown in the other figures in this section shows that when 
parameters in the waveform model other than M and R must be determined as part of the waveform 
fitting process, the uncertainties in the resulting estimates of M and R are generally much larger than 
when only M and R are unknown. The reason the uncertainties are so large is that the waveform is 
degenerate with respect to changes in many of the parameters in the waveform model, i.e., the effects 
on the waveform of changing the value of one of the system parameters can be partially or wholly 
compensated by changing the value of one or more other parameters: 

— The effect on the stellar compactness and surface redshift of an increase in R can be compensated 
almost exactly by an increase in M, because the special relativistic effects on the waveform of the 
rotational surface velocity are ~ 5% or less for the stellar models considered here (see Section [3j). 
This is a very strong degeneracy. 

— Changes in the spot inclination, viewer inclination, and stellar radius affect the waveform in similar 
ways and can therefore compensate for one another. The increase in the relativistic Doppler boost 
and aberration produced by an increase in the spot inclination can be compensated in part by a 
decrease in the observer's inclination and/or a decrease in the stellar radius. The change in the 
projected emitting area seen by the observer produced by an increase in the spot inclination can 
be partially compensated by an increase in the observer's inclination. 

— An increase in the stella r compactness can b e partially compensated for by an increase in the spot 



radius, but as shown by lLamb et al.l (|2009al ). the dependence of the waveform on the spot radius 
is weak unless the spot is very large, while its dependence on the compactness is weak unless the 
star is very compact {M/R > 0.25). Hence this is a weak degeneracy. 

Changing the number and spectrum of the background counts affects the inferred spectrum and 
modulation fraction of the waveform and therefore can partially compensate for changes in any of 
the other parameters that affect the spectrum and modulation fraction of the waveform, including 
the spot inclination and size, the stellar compactness, and the observer inclination. 
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As noted previously, these degeneracies are an inherent property of any physical model based on a 
rotating hot spot and cannot be removed by "improving" the model. 

As a result of these degeneracies, the uncertainties in estimates of M and R depend sensitively on 
the physical properties of the system. If the hot spot and the observer are at high inclinations, the 
effects on the waveform caused by the rotation of the star will be maximal, reducing the effects of these 
degeneracies. If instead the hot spot and the observer are at low inclinations, the effects caused by the 
rotation of the star will be much less and these degeneracies will much more important. 

A large background exacerbates the degeneracy problem, because it increases the statistical fluctu- 
ations in the observed waveform. As a result, a wider range of models will adequately fit the waveform 
data. On the other hand, independent additional knowledge about the system can reduce or eliminate 
degeneracies. For example, accurate a priori knowledge of the observer's inclination can significantly 
improve the constraints, if the spot and observer inclinations are high; a priori knowledge of the dis- 
tance to the system can also help. Independent knowledge of the background can be used to tighten the 
constraints on M and R substantially. If the spot and observer inclinations are high, measurement of an 
identifiable atomic spectral line in the burst oscillation spectrum can be used to tighten the constraints. 



4-2.3. Effects of spot and observer inclinations 

Figure E] (see also Table [2]) shows that the constraints on M and R are fairly tight when the hot spot 
and the observer are both at high inclinations but very loose when they are both at low inclinations, 
regardless of the background level. Comparison of Figures Eh, Eti, and [3b (again see Table [2]) shows 
that for high observer inclinations, the constraints on M and R weaken slowly as the spot inclination 
decreases from 90° to 60°. 

If the hot spot and observer inclinations are both high, the la uncertainties in M and R range from 
a few percent for our low and medium backgrounds to about 10% for our high background (see Table [2]). 
The 3(T uncertainties are much larger, ranging from about 15% for our low or medium backgrounds to 
about 50% for our high background. If instead the spot and observer inclinations are both low, the 
Icr uncertainties in M and R are large, ranging from 30% to 50% or more (again see Table [2]), and 
the 3o" contours span much or all of the M-R domain that was searched. The importance of high 
spot and observer inclinations becomes especially clear if one compares the confidence regions for the 
high- inclination case shown in Figure [2]i with the confidence regions for the low- inclination case shown 
in Figure [2^. The fractional rms amplitudes of the oscillations are the same in these two cases; only the 
inclinations are different. 

For high hot spot and observer inclinations, the oscillation amplitude and the effects on the energy 
dependence and harmonic content of the waveform of the relativistic Doppler boost and aberration are 
responsible for the lower bound on R. The upper bound on R is provided by the requirement that the 
amplitude be sufficiently large but the effects of rotational motion on the waveform must not be too 
large. This requirement also imposes a lower bound on the observer's inclination. 

For low spot and observer inclinations, the waveform is nearly sinusoidal and the modulation 
fraction is low. Therefore much wider ranges of the parameters other than M and R (such as the spot 
and observer inclinations) give a modulation fraction low enough to be consistent with the observed 
waveform, if the compactness is high. Consequently, when we marginalize over these other parameters. 
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M-R pairs near the high compactness boundary are given more weight. As a result, the highest posterior 
probabihty density is typically near the highest compactness that was searched, which means it is 
typically close to the high-compactness boundary of the search domain, even though comparable fits 
would be possible for many values of M/R if we were to optimize the other parameters. 



4-2.4- Effects of the background 

Figure [2] (see also Table [2|) illustrates how the constraints on M and R are affected by background 
counts. The background has very little effect on the uncertainties unless the number of background 
counts exceeds the number of counts from the hot spot, in which case the uncertainties increase as 
the background increases, approximately as 1/7^, where IZ = A^osc/V-^tot in terms of the number of 
counts Nose in the oscillation and the total number of counts Atot (see equation ([1|)). Thus, the la and 
3cr uncertainties for our high background, which contributes 9 times as many counts as the hot spot, 
are several times larger than the uncertainties for our low and medium backgrounds, which contribute, 
respectively, less than and about the same number of counts as the hot spot. 

The values of IZ and the fractional rms amplitudes of the oscillations in the synthetic waveforms 
generated in the low- and high-inclination reference cases for our low, medium, and high backgrounds are 
listed in Table [H Scaling from RXTE observations suggests IZ ~ 150 for a typical burst oscillation ob- 
served during a single burst using a 10 m^ detector (see Section [2.2. ip . TZ scales as the square root of the 
number of counts, so combining the data from 25 such bursts (see Section [4.2.71 and Appendix[B]) could 
give 7Z ~ 800, between the values for the synthetic waveforms generated in the high-inclination, medium- 
background and high-inclination, high-background cases. The high-inclination, high-background case 
has a fractional rms amplitude comparable to the ~ 10% fractional amplitudes typically observed in 
burst tails using RXTE. 



4-2.5. Effects of independent knowledge of relevant system properties 

The extent to which knowledge that is independent of the waveform fitting process improves the 
precision of M and R estimates depends strongly on the inclinations of the hot spot and the observer. If 
they are high, the uncertainties in M and R can be reduced by independent knowledge of the observer's 
inclination or of the distance to the star. In this case observing an atomic scattering line with a known 
rest energy tightly constrains the stellar compactness M/R and improves the constraints on M and R 
more than knowing the distance and approximately as much as knowing the observer's inclination. If 
instead the spot and observer inclinations are low, observing a scattering line with a known rest energy 
still tightly constrains the stellar compactness but improves the constraints on M and R very little. 
The uncertainties in M and R can be reduced substantially by independent knowledge of the size and 
spectrum of the background, especially if the spot and observer inclinations are high. Table [2] lists the 
Icr fractional uncertainties for each of these cases. Analysis of high-quality observations of the spectra 
of the bursts using high-precision spectral models can help to further constrain M and R. 

Figured^ shows that if the spot and observer inclinations are both 90°, knowing the inclination 
of the observer reduces the Icr uncertainty from ~ 9% to ~ 5% and the 3cr uncertainty from ~ 50% to 
~ 20% for our high background (compare this figure with Figure [2^) . Figure Hb shows that knowing 
the distance to the star reduces the Icr uncertainty to ~ 7% and the 3cr uncertainty to ~ 30% (again 
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Table 3. Scaling of uncertainties in M and R with total counts'' 



iVtot 


5Mi{%) 


SMs{%) 


dRi{%) 


SR3{%) 


10^ 


15 


42 


12 


27 


10^ 


6.8 


15 


4.7 


11 


106 


1.8 


5.0 


1.3 


3.1 



a/ 



'A'^tot is the total number of counts in the high- 
inclination, medium-background synthetic waveforms that 
were analyzed, 5Mi and 5M-^ are the la and 3a uncertain- 
ties in M, and 5Ri and 6R3 are the la and 3a uncertainties 
in R. These uncertainties assume all system parameters 
are known except M and R (see text). The fractional rms 
amplitudes of all three waveforms are ~ 54% (see Table [3|. 



Table 4. TZ values and rms amplitudes for selected synthetic waveforms'' 



Inclination 


Background 


n 


rms 


low 


low 


340 


0.21 


low 


medium 


280 


0.14 


low 


high 


120 


0.028 


high 


low 


1360 


0.81 


high 


medium 


1100 


0.54 


high 


high 


490 


0.11 


''Here 7^ = 


A'osc/\/A'tot, where iVosc is the 



number of counts in the oscillating compo- 
nent of the synthetic waveform and A'tot is 
the total number of counts in the waveform. 
See equation ([T]). 
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compare this figure with Figure [2]3). Additional computations show that if the spot is at 80° and the 
observer's inchnation is known, the uncertainties in M and R are similar to the uncertainties obtained 
if the spot is at 90°. If the spot and observer inclinations are high and the observer inclination and the 
distance are both known, the la uncertainties in M and R are similar to the uncertainties obtained if 
only the observer's inclination is known. The similarity of the la uncertainties in these two cases is 
not surprising, because independent knowledge of the observer's inclination is usually more helpful than 
independent knowledge of the distance. If instead the spot and observer inclinations are low, knowing 
the inclination of the observer or the distance to the star reduces the uncertainties only modestly 
(compare Figures Sja and [Hi with Figure [Sf) . 

Figures 0^ and Hf illustrate the effect on the precision of M and R estimates if an atomic scattering 
line with a known rest energy is observed in the emission from the hot spot. The results shown here 
are for a line at 6.4 keV with a FWHM of 0.2 keV and an optical depth at the line center of 0.24. If the 
spot and observer inclinations are high, observing this line tightly constrains the stellar compactness 
and tightens the constraints on M and R, decreasing the la uncertainty from ~ 9% to ~ 4% and the 
3a uncertainty from ~ 50% to ~ 20% for our high background, a greater improvement than would be 
achieved by independently determining the distance and approximately as much as would be produced 
by independently determining the observer's inclination (compare Figure 0^ with Figure [2^). If instead 
the spot and observer inclinations are low, observing this scattering line still tightly constrains the stellar 
compactness but improves the constraints on M and R very little (compare Figure Hf with Figure [2]f; for 
the particular realization shown in Figure Sf, the constraints on M and R are weaker with a scattering 
line present, probably due to a sampling fluctuation). 

Figure [5] shows that the knowing the size and spectrum of the background greatly improves the 
constraints on M and R. For our high-inclination case, knowing the background decreases the la un- 
certainties in M and R from ~ 9% to ~ 4% and the 3a uncertainties from ~ 50% to ~ 9% (compare 
Figure [5^ with Figure [2^). For our low- inclination case, knowing the background leads to a la un- 
certainty of ~20%, whereas without this knowledge one obtains no useful constraints on M and R 
(compare Figure \5}p with Figure [2]F) . 

Fitting high-quality observations of the X-ray spectra of many bursts from the same star using 
detailed, high-precision model atmosphere calculations can produce a constraint on M and R that is 
complementary to that obtained by fitting burst oscillation waveforms. This is illustrated by the dash- 

2 

dotted red curves in Figures[2^,[3^,[3)3,[ll;,[l^,ll]F,[6^, and[6)3, which show the relation {l + z)/g9 = const. 
for the values of the redshift z and surface gravity g that correspond to the M and R values used in 
generating these synthetic waveforms. As discussed in Section l4.1.3l this combination of g and z can be 
determined by accurately measuring the spectra of the X-ray bursts. The red dashed-dot curve intersects 
the confidence region derived by fitting burst oscillation waveforms at a high angle. Combining these 
two complementary methods is therefore a promising way to further constrain M and R. 



4-2.6. Effects of errors in the assumed properties of the hot spot 

In the results discussed previously, we assumed that we knew the shape of the hot spot and the 
spectrum and beaming pattern of the radiation from the spot. In practice, we may not know these 
properties. Any discrepancies between the actual properties of the hot spot and the properties assumed 
in the waveform model will produce systematic errors in the constraints on M and R, in addition to 
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the statistical uncertainties, and these errors should be quantified. A full analysis of possible systematic 
errors is beyond the scope of this work, but our initial exploration shows that inaccuracies in modeling 
the waveform can sometimes increase the uncertainties in measurements of M and R. 

Errors in the assumed spot shape. The effects of errors in the assumed shape of the hot spot are 
illustrated in Figured For these examples, we generated synthetic waveforms using hot spots elongated 
in the latitudinal (east-west) and longitudinal (north-south) directions by superimposing two or three 
circular emitting areas, as described in the figure caption, and then fitting the resulting waveforms using 
our model, which assumes a circular hot spot. 

Figure [6^ shows the results for the moderately high-inclination spot elongated in the north-south 
direction by 50° and our medium background. The la and 3cr uncertainties in M and R are, respectively, 
~40% and ~75%, much larger than the ~4% and ~ 15% uncertainties for the waveform produced by 
a circular spot in the rotational equator and the same background (see Figure ^). The fit appears 
excellent (x^/dof = 432.5/443), providing no warning that the model is incorrect. Although the best- 
fit values of M and R differ substantially from the values assumed in generating the waveform, the 
differences are only marginally significant. 

Figure [lb shows the results for the spot in the rotational equator elongated in the east-west direction 
by 45° and our medium background. The Icr uncertainties in M and R are, respectively, ~ 8% and 
~ 17%, substantially larger than the ~4% uncertainties for the waveform produced by a circular spot 
in the rotational equator and the same background (see Figure [2}:) . The 3a uncertainties in M and R 
are, respectively, ~30% and ~40%, much larger than the ~ 15% uncertainties that were achieved for 
a circular spot (again see Figure [2}:). The fit is good (x^/dof = 442/443), providing no warning that 
the model is incorrect. However, in this case the best-fit values of M and R are very close to the values 
assumed in generating the waveform. 

Errors in the assumed spot spectrum. The effects of errors in the assumed spectrum of the emission 
from the hot spot are illustrated in Figures [7^. and Eh- For these examples, we generated synthetic 
waveforms assuming the spectrum of the radiation from the surface of the hot spot has the shape 
of a Bose-Einstein spectrum and then fitted the waveforms using our model, which assumes that the 
radiation has the shape of a Planck spectrum. 

For the high-inclination reference case (Figure [7ti) , the best-fit values of M and R are very close 
to their actual values and the quality of the fit is acceptable (there would be a 10% probability of a 
reduced x^ this high if the model were correct). The la uncertainties in M and R are both about 6%, 
whereas the 3a uncertainties in M and R are 15% and 12%, respectively. These are comparable to the 
3%-4% la and about 15% 3a uncertainties achieved by fitting a model with the correct spectral shape 
to the synthetic waveform (see Section [4. 2. 3p . 

For the low-inclination reference case (Figure [7b), the confidence regions are large, the la region 
extends beyond the upper boundary of the search domain, and the best-fit M-R pair is on the boundary, 
all indicating that the most probable solution has not been determined. The reduced x^ value is close to 
unity and hence does not itself show that the waveform model is wrong, but this is of little consequence 
because the constraints on M and R are very weak. 

Errors in the assumed spot beaming function. The effects of errors in the assumed beaming pattern 
of the radiation from the hot spot are illustrated in Figures [7h and[7ll. For these examples, we generated 
synthetic waveforms assuming that the surface of the hot spot radiates isotropically, but then fit these 
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waveforms using our standard model, which assumes the beaming pattern from the surface of the hot 
spot is described by the Hopf function (the beaming pattern for radiation from a Thomson scattering 
atmosphere) . 

For the high-inchnation reference case (see Figure Eb), the la and 3a contours are very small and 
far from the M-R pair used in generating the synthetic waveform, but the reduced x^ for this fit is 
extremely large, showing that the model being fit is wrong. 

For the low-inclination reference case (see Figure [TJi) , the best-fit M-R pair is again far from the 
pair used in generating the synthetic waveform. Although the reduced x^ value is close to unity and 
hence does not itself show that the waveform model is wrong, this is of little consequence, because the 
constraints on M and R are very weak, the la and 3a confidence regions extend beyond one or both of 
the boundaries of the search domain, and the best-fit M-R pair is on the high-compactness boundary, 
all indicating that the most probable solution has not been determined. 



4-2.7. Combining different segments of data 

We have seen that for favorable system properties, obtaining strong constraints on M and R usually 
requires ~ 10^ counts from the hot spot. Acquiring this many counts will typically require collection 
of hundreds to thousands of seconds of data, even using an instrument with an effective area ~ 10 m^. 
Thus, unless a superburst is observed, the required data must be accumulated over many bursts. We 
therefore analyzed how the constraints on M and R reported here would be affected if ~ 10® hot spot 
counts were collected by combining data from multiple bursts from the same neutron star. The same 
considerations apply if the data from a single burst are divided into several time segments that are then 
analyzed jointly. 

We find that the constraints on interesting parameters obtained by jointly fitting many data sets 
are usually comparable to the constraints that could be obtained by fitting a single set of waveform 
data with the same average profile and the same total number of counts as the waveform obtained by 
combining the data sets (see Appendix[B]). Thus, if '^ 10® counts can be accumulated from the hot spots 
of many bursts produced by a given star, the prospects are good that constraints on M and R similar 
to those reported here can be obtained. This is true even though one might intuitively expect that the 
extra parameters involved in fitting many data sets with changing waveforms would compromise the 
constraints on the model parameters that are fixed (M, R, and ^obs)- 

The effects on the waveform of the special relativistic Doppler boost and aberration play an im- 
portant role in constraining M and R. These effects can be accurately captured even if the hot spot 
rotation frequency is changing, provided it can be tracked accurately enough to maintain the correct 
oscillation phase when folding successive periods of the oscillation. If this can be done, the constraints 
on M and R that can be obtained by analyzing the resulting folded waveform are often nearly the same 
as those that could be obtained by analyzing a similar waveform with a fixed oscillation frequency and 
the same number of counts; moderate changes in the values of nuisance parameters other than the spot 
rotation frequency often have only a small effect. 

If the oscillation frequency varies so rapidly or erratically during a burst that it cannot be tracked 
or other important parameters (such as the spot inclination and spot radius) vary sufficiently during 
a single burst or from burst to burst that this has a substantial effect on the waveform, the full data 
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set may have to be divided into a series of shorter segments and each segment analyzed separately. We 
find that the resulting constraints on M and R will usually be similar to the constraints obtained by 
analyzing a single data set in which these nuisance parameters do not vary and that has the same number 
of counts. Analyzing a sequence of data segments increases the computational burden only linearly with 
the number of segments. If the segments are properly analyzed using a Bayesian approach, the order 
in which they are analyzed does not matter. 

In summary, constraints similar to those reported here can usually be obtained if 10^ hot spot 
counts are acquired by combining data from different segments of a single burst or from multiple bursts 
from the same neutron star, even if the oscillation frequency and other physical parameters vary during 
the burst or from burst to burst. 
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Fig. 2. — These results show that M and R are tightly constrained when the spot and observer incli- 
nations are high and the total number of counts from the background is comparable to or less than the 
number from the hot spot. The contours in the left column are for the high spot and observer inclination 
model and (top to bottom) the low-, medium-, and high-background models (for the properties of these 
models, see Tables [1] and [2|) . The right column shows the corresponding contours for the low spot and 
observer inclination model. The dashed lines show the R/M = 3.5 and R/M = 6.8 boundaries within 
which the posterior probability distribution was computed. The dotted and solid curves show, respec- 
tively, the portions of the la and 3a confidence contours within this domain. The + symbol indicates 
the radius and mass that were used in generating the synthetic waveform and the star symbol indicates 
the values where the marginalized posterior probability density is highest. The red dashed-dot curve 
in Fig. [2^ shows points of constant (79/(1 + z), illustrating the complementary constraint that can be 
obtained by fitting the spectra of X-ray bursts from the same star (see text). 
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Fig. 3. — Comparison of the contours shown here and in Fig. [2}: demonstrates that the constraints on 
M and R weaken slowly as the spot inclination decreases from 90° (in Fig. [2]3) to 80° (left panel) and 
to 60° (right panel). All three plots assume the observer inclination is 90° and our medium background 
(~ 10^ counts from the hot spot and ~ 10^ background counts). The red dashed-dot curves show points 

2 

of constant 9 9 /(I + z\ illustrating the complementary constraint that can be obtained by fitting the 
spectra of bursts from the same star (see text). For the meanings of the line types and symbols, see 
Fig. El 
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Fig. 4. — These results show that independent knowledge of some system parameters can significantly 
tighten the constraints on M and R. All results are for our high background. The contours in the 
left column are for the high spot and observer inclination case and (top to bottom) known observer 
inclination, known distance, and an identifiable atomic scattering line in the waveform spectrum (see 
text for details); these contours should be compared with those in Fig. [2^. The right column shows 
the corresponding contours for the low spot and observer inclination case; these contours should be 
compared with those in Fig. [2f. For the meanings of the line types and symbols, see Fig. [2j 
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Fig. 5. — These results assume that the properties of the background are known. Left: Constraints for 
the high-inchnation reference case (spot and observer inchnations of 90°). Right: Constraints for the 
low-inchnation reference case (spot and observer inchnations of 60° and 20°, respectively). Comparison 
of these contours with those shown in Figs. ^ and [2li, which assume the background is unknown, 
demonstrate that M and R are much more tightly constrained if the background is known. All these 
results assume our medium background. For the meanings of the line types and symbols, see Fig. [2j 
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Fig. 6. — These results show that fitting the observed waveform using a model with an incorrect spot 
shape increases the uncertainty of the M and R estimates and can sometimes bias them. Left: Con- 
straints obtained when the actual emitting area is elongated in the north-south direction, represented 
here by two uniform, circular spots of angular radius 25°, centered at the same longitude, with incli- 
nations of 40° and 90°, but the waveform is fit assuming a circular spot. Right: Constraints obtained 
when the actual emitting area is rotational equator and elongated in the east-west direction, represented 
here by three uniform, circular spots of angular radius 25° , with centers at the same 90° inclination and 
spaced 22.5° apart in the azimuthal direction, but the waveform is fit assuming a circular spot. Both 
results are for our medium background. For the meanings of the line types and symbols, see Fig. [2j 
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Fig. 7. — The top panels show that fitting the observed waveform using a model in which the emission 
spectrum differs by a modest amount from the actual emission spectrum has little effect on the estimated 
values of M and R or their uncertainties (compare Figs. [2^ and [2)3). The bottom panels show that 
fitting the observed waveform using a model in which the beaming pattern differs substantially from the 
actual beaming pattern increases the uncertainty of M and R estimates and can sometimes bias them a 
moderate amount. The contours in the left column are for high spot and observer inclinations, whereas 
the contours in the right column are for low spot and observer inclinations. All results assume our 
medium background. The top row shows the constraints obtained when the actual emission spectrum 
has the shape of a Bose-Einstein spectrum with a chemical potential /i = —kT but the waveform 
model assumes the emission spectrum has the shape of a Planck spectrum. The bottom row shows 
the constraints obtained when the actual emission is isotropic but the waveform model assumes the 
emission has the beaming pattern described by the Hopf function. For the meanings of the line types 
and symbols, see Fig. [2l 
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5. SUMMARY AND CONCLUSIONS 

The goal of this paper was to explore the constraints on the masses and radii of neutron stars that 
could be inferred by analyzing observations of X-ray burst oscillations made using a future, satellite- 
borne detector with 2-30 keV energy coverage and an effective area 10 to 20 times larger than the RXTE 
PC A. This research was motivated by the new concepts for large-a, r ea X-ray timing space missions that 



are now being proposed. These include LO FT (JMignani et al.ll2012l : lDel Monte et al.ll2012l ) and AXTAR 



(jChakrabartv et alj 120081 : iRav et alj 120111 ) . Although we focused here on exploring the constraints on 
M and R that can be obtained by analyzing the waveforms of X-ray burst oscillations, our results are 
equally relevant to plans for analyzing the oscillations produced by X-ray emission from the heated 
polar caps of isolated rotation-powered millisecond pulsars, the goal of the proposed NICER mission 



(JGendreau et alJl2012l ). 



Approach to the problem. We described and explained our approach in Section [2j There we 
considered the relative merits of using observations of oscillations during burst rises and burst tails (see 
Section [2. 2p . We find that the uncertainties in M and R estimates obtained by analyzing oscillations 
during the tails of bursts are likely to be smaller than the uncertainties in estimates derived by analyzing 
the oscillations observed during the rises of bursts. Hence analyzing tail oscillations is likely to be the 
better approach. 

Computational methods. We described our computational methods in Section [3l We explored the 
constraints on M and R that could be derived from X-ray burst oscillations by first generating energy- 
dependent synthetic observed waveforms for a variety of neutron star and hot spot properties. We then 
used a Bayesian approach and MCMC sampling methods to determine the constraints on M and R that 
can be obtained by analyzing these synthetic observed waveforms. We computed the joint posterior 
probability distribution of the parameters in our waveform model, marginalized the parameters other 
than M and R, and then determined the most probable values of M and R and their confidence intervals. 

The purpose of the analysis presented here was not to reproduce all the steps that would be needed 
for a full Bayesian analysis of a real observed waveform, but rather to determine the precision and 
accuracy with which such a full analysis could determine M and R. We therefore made use of several 
shortcuts that reduced the computational burden substantially without significantly altering the results. 
Using these shortcuts, we could compute the most probable values of M and R and their confidence 
intervals for a single synthetic observed waveform in 50-100 clock hours, running the parallel version of 
our code on 150 nodes of a fast CPU cluster. 

Assum,ed counts from the hot spot and background. We assumed that about 10^ hot spot counts 
are available from the star. A 10 m^ detector could collect this many hot spot counts by observing 
20-25 bright X-ray bursts (see Section 12.20 . We find that the constraints on M and R obtained by 
combining data from multiple segments of a single burst or from multiple bursts from the same star 
are usually similar to the constraints that would be obtained by analyzing an observation of a single, 
time-independent waveform that has the same number of counts and the same shape as the average of 
the waveform over the multiple data sets (see below, and Section [4.2.71 and Appendix [B]) . 

We assumed that all sources of background, when combined, contribute about 0.3 x 10^, 10^, or 
9 X 10^ counts. Our treatment of the background was very conservative, in the sense that we usually 
made no assumptions about the magnitude or spectrum of the background. We did not even assume 
that the background is constant, only that it does not vary at frequencies commensurate with the burst 
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oscillation frequency. 

Our main results. We described our results in detail in Section [H There we showed that if ~ 10^ 
counts were collected from the hot spot of a single bursting neutron star and the values of all the 
parameters in the waveform model other than M and R were known independently of the waveform 
analysis, the uncertainties in M and R estimates produced by statistical fluctuations in t he waveform 



would typically be only a percent or two (compare the analysis of a similar situation in IStrohmaver 



20041 ) . For realistic situations in which some of the other parameters in the model must be determined 
from the waveform, the statistical uncertainties are usually much larger. The reason for this is that the 
effects on the waveform of changing different parameters in the waveform model are often very similar. 
These degeneracies of the waveform with respect to changes in the values of the model parameters are 
an inherent property of any physical model based on a rotating hot spot and cannot be removed by 
"improving" the model. These degeneracies inflate the uncertainties in M and R estimates produced 
by the statistical fluctuations in the waveform. 

We first explored how the uncertainties in the measured values of M and R depend on the inclina- 
tions of the hot spot and observer and the background count rate. We find that the uncertainties depend 
strongly on the inclination of the hot spot and observer relative to the stellar spin axis. Specifically: 

— If the hot spot is within 10° of the rotation equator, both M and R can usually be determined 
with an uncertainty of about 10%. 

— If instead the spot is within 20° of the rotation pole, the uncertainties in M and R are so large 
that waveform measurements alone provide no useful constraints. 

— The uncertainties in M and R are affected little by background count rates less than or comparable 
to the count rate from the hot spot, but become significantly larger for higher background count 
rates. 

Next, we explored the effects on the uncertainties in M and R if the distance to the star or the 
inclination of the observer are known from other measurements, if a resonance scattering line is observed 
in the burst oscillation spectrum, or if the properties of the background are independently known. We 
find that: 

— Independent information about the background can greatly reduce the uncertainties in estimates 
of M and R. 

— Independent knowledge of the observer's inclination can also greatly reduce the uncertainties. 

— Knowledge of the star's distance can help, but not as much as knowledge of the background or 
the observer's inclination. 

— Observation of an identifiable atomic line in the hot-spot emission always tightly constrains M/R; 
it can also tightly constrain M and R individually, if the hot spot is within about 30° of the 
rotation equator. 

Finally, we explored the effects of deviations in the actual shape of the hot spot, radiation beaming 
pattern, and spectrum from those assumed in the our waveform model. We find that: 

— Modest deviations of the actual spectrum from that assumed in the waveform model have little 
effect on the accuracy or uncertainty of M and R estimates. 
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Large deviations of the actual radiation beaming pattern from the pattern assumed in the wave- 
form model can increase the uncertainties of M and R estimates substantially. 

In some cases, but not always, large deviations of the actual shape of the hot spot from the circular 
shape assumed in the waveform model can increase the uncertainties of M and R estimates and 
bias them by moderate amounts. The physical conditions that produce tight constraints on M 
and R (relatively small spots far from the rotat ion axis) are the conditions in which the shape of 
the spot is unimportant (see lLamb et al.l l2009a| ]bl) . 



Combining data from different bursts. Our results show that for favorable system properties, strong 
constraints on M and R can be obtained if one can analyze ~ 10^ counts from the hot spot. This will 
typically require analyzing hundreds or thousands of seconds of data, even if the detector has an effective 
area ~ 10 m^. Thus, unless a superburst is observed, the data required must be accumulated from many 
bursts. We therefore analyzed how the constraints on M and R reported here would be affected if the 
~ 10^ hot spot counts were obtained by combining data from multiple bursts from the same star. The 
same considerations apply if the data from a single burst are divided into several time segments and 
then analyzed jointly. We find that the constraints on M and R obtained by jointly fitting many data 
sets are usually comparable to the constraints that could be obtained by fitting a single segment of 
waveform data having the same profile and as the average profile of the combined data sets and the 
same total number of counts (see Section 14.2.71 and Appendix [B]) . 

Time varying waveforms. We also investigated the problem of analyzing burst oscillation waveforms 
that change with time. If the oscillation frequency changes during a burst but the change can be 
modeled accurately enough to maintain the correct oscillation phase when folding successive periods 
of the oscillation, the constraints on M and R obtained by analyzing the resulting folded waveform 
will be nearly the same as those that could be obtained by analyzing a similar waveform with a fixed 
oscillation frequency and the same number of counts. Even if the oscillation frequency varies too rapidly 
or irregularly during the burst rise or tail to be described accurately by a simple frequency model, the 
full burst oscillation data set can be divided into smaller time segments and analyzed using standard 
Bayesian techniques. This approach can also be used if other physical properties of the system, such as 
the size and inclination of the emitting region, vary significantly. The computational burden of this kind 
of analysis increases only linearly with the number of segments. Consequently, variations of the burst 
oscillation waveform on timescales shorter than the burst rise or burst tail (but substantially longer 
than the burst oscillation period) do not appear to pose an insurmountable analysis problem (again see 
Section [3221 and Appendix |B]) . 
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CODE TESTS 



In this appendix, we summarize the suite of test problems we routinely solve to validate both our 
waveform code and our MCMC code (of which the waveform code is an essenti a l part) . The waveform 



code we use here is based on the code developed and validated bv iLamb et al.1 (|2009al ). We described 
some tests of this code there. The code validation tests we describe here are designed to check many of 
the physical effects included in our codes, by using them to solve problems where the answer is known 
analytically. This suite of test problems includes the following problems: 

— computing the deflection angle as a function of the star's compactness and the direction of the 
ray relative to the normal to the stellar surface (Section lA.l.ip 

— computing the light travel time as a function of the star's compactness and the direction of the 
ray relative to the normal to the stellar surface ()A.1.2p 

— computing the bolometric flux seen by a distant observer located directly above the center of an 
emitting spot on a non-rotating star (Section IA.1.3P 

— computing the radiation pattern produced by a small emitting spot on a non-rotating star in flat 
spacetime (Section lA.1.4p 

— computing the waveform produced by a small emitting spot on a rotating star in flat spacetime 
(Section lA.l.Sp 

— computing the angular deflection of a pencil beam from a small emitting spot on a rotating star 
(Section IAXG]) 

— computing the observed temperature of thermal emission from a small emitting spot on a rotating 
star in flat spacetime, as a function of rotational phase (Section IA.1.7P 

— determining the stellar compactness, spot temperature, and distance, given uniform thermal emis- 
sion from the entire surface of a non-rotating star with no background (Section IA.2.ip 

— determining the spot inclination, observer inclination, and distance, given thermal emission from 
a small hot spot on a non-rotating star in fiat spacetime with no background (Section IA.2.2P 

— determining the stellar compactness, spot temperature, and distance, given thermal emission from 
a small spot on a non-rotating star in flat spacetime with a background (Section lA. 2. 3p . 

The values of the waveform code integration, hot spot, and angular resolution parameters that we 
use for the code tests reported here are the values listed in Table lAll unless otherwise noted. We use 
these same parameter values when computing M and R confidence regions (see Section 14. 2p , except 
that for that purpose we set Nphase to 16, which speeds up the computations while providing accuracy 
sufficient for the waveform-fitting process. As discussed in Section IA.31 we have carefully chosen the 
values of these code parameters to provide a resolution fine enough to meet our accuracy requirements, 
but no finer, so that our code runs as fast as possible. 

We now discuss the suite of validation tests in detail. 
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A.l. Tests of the Waveform Code 



A. 1.1. Deflection angle as a function of compactness and direction 

In the Schwarzschild spacetime, the total angle il^ by which a light ray from the stellar surface is 
eventually deflected is a function only of the initial angle a of the ray relative to the normal to the 
stellar surface (as i neasured in the static frame) and the compactness M/R of the star, and is given by 
the expression (see iPechenick et al.lll983l . equation (2.13) for a similar expression) 



/ / ^ 



sin a dx 



\/(l - 2M/R) - (1 - 2Mx/R) x^ sin^ a 



(Al) 



The integrand in equation (lAlh diverges as x — )• 1 and sin a — )• 1; numerical integration near this limit 
therefore requires special treatment. In the waveform code we break the integral into two pieces, from 
x = Otox = l — e^ and from x = 1 — ed to x = 1^ with ej^ <C 1. Our waveform code computes the 
first piece numerically, using Simpson's rule with n^ divisions, and computes the second piece using an 
approximate analytical expression, derived as follows. 

Setting X = \ — e^, the square root in the integrand becomes 
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Dropping the 0{e^ terms inside the square root, the second piece of the integral becomes 
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where a, 6, and c are the coefficients of 1, ed, and e^ in the square root (jA3J). Because c = sin^ a {6M/R — 1) 
different integrals are performed for R < 6M and R > 6M. Note also that for R = 3Af and sin a — )• 1, 
a = 6 = 0, causing the integral to diverge, as it must, because r = 3M is the photon orbit in the 
Schwarzschild geometry. 

We test the deffection angle '0 given by our waveform code by comparing it with the value of 
tp given by evaluating integral (lAlh using Mathematica 8. Figure lAll shows the fractional difference 
[{tpcode - '0mathematica)/'0mathematica] for three valucs of Ud and two valucs of the Stellar compactness. 
Figure [A2] shows the fractional difference for three values of Cd and two values of the stellar compactness. 
The left and right panels of Figure lA3l show, respectively, the relation ^l'{a) and the fractional difference 
in ■0 given by our code using Ud = 100 and €d = 0.1, for four values of the stellar compactness that span 
the range considered in this work. Figures IA1HA3I show that the values Ud = 100 and e^ = 0.1 of these 
resolution parameters that we use in our waveform code when fitting synthetic waveform data provide 
sufficient accuracy for this purpose. 
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Fig. Al. — Fractional difference between the deflection angle ^ as a function of the initial ray direction 
a. computed using our code and using Mathematica 8, for three different numbers n^ of integration steps 
and two values of the stellar compactness (left: R/M = 5.0; right: R/M = 4.0). This figure shows that 
the value n^ = 100 used in our waveform fitting is adequate. See Section [A. 1 . 1 1 for further details. 
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Fig. A2. — Fractional difference between the deflection angle i/j as a function of the initial ray direction 
a computed using our code and using Mathematica 8, for three values of the integration parameter e^ 
and two values of the stellar compactness (left: R/M = 5.0; right: R/M = 4.0). This figure shows that 
the value e^ = 0.1 used in our waveform fitting is adequate. See Section [A. 1.1 1 for further details. 

A. 1.2. Light travel time as a function of compactness and direction 

The time it takes for a light ray emitted at an angle a relative to the surface normal (in the static 
frame) to travel from the surface to infinity, as measured at infinity, is given by 
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This expression accounts for the gravitational time delay in the Schwarzschild spacetime as well as the 
purely geometrical effect that the travel time is different for rays emitted from the stellar surface at 
different angles. The travel time given by equation ()A5P is formally infinite, because the ray is traced 
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Fig. A3. — Lejt: Deflection angle ^ as a function of the initial ray direction a. computed using our code, 
for four values of the stellar compactness. The deflection angle is larger for smaller values of R/M, i.e., 
for more compact stars. Right: Fractional difference between the deflection angle ip computed using our 
code and using Mathematica 8, as a function of a, for the same four values of the stellar compactness. 
The fractional difference is larger for smaller values of R/M, i.e., for more compact stars, but is smaller 
in magnitude than 2.5 x 10~^ for all these values of R/M. This figure shows that the resolution 
parameter values rid = 100 and e^ = 0.1 used here provide sufficient accuracy for our waveform fitting. 
See Section lA.l.ll for further details. 



to infinity. However, what we are interested in is not the total travel time of each ray, but rather the 
differential travel times of different rays. 

To be able to compute the differential travel times of any two rays, we calculate the differential 
travel time of every ray relative to the travel time of a reference ray. We choose as our reference ray 
the radial ray coming from the point on the stellar surface directly under the observer. The differential 
light travel time of an arbitrary ray compared to the reference ray is then 
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The integrand in this expression can diverge at the two limits of the integral. Our waveform code 
therefore calculates this integral in three pieces. Near the upper endpoint, we set y = 1 — ei, with 
ei <C 1, expand the expression inside the square root in powers of ei, and drop the terms that are 0(ef ). 
The result is 
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where z = 1 - 2ei + ei(2M/i?)/(l - 2M/R) and zi = 1 - 2ei(l - 3M/R)/{1 - 2M/R). Near the lower 
endpoint of the integral, we set y = €q, with eo ^ 1, expand the integrand in powers of eo, and thereby 
obtain the expression 
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Fig. A4. — Fractional difference between the differential light travel time At as a function of the initial 
ray direction a computed using our code and using Mathematica 8, for three different numbers nt of 
integration steps and two values of the stellar compactness (left: R/M = 5.0; right: R/M = 4.0). This 
figure shows that the resolution nt = 100 used in our waveform fitting is adequate. See Section IA.1.21 
for further details. 



Our waveform code computes the middle piece of the integral, from y 
using Simpson's rule with nt divisions in this interval. 



eo to y = 1 — ei, numerically, 



We test the differential light travel time At given by our waveform code by comparing it with the 
value of At given by evaluating integral ()A6P using Mathematica 8. Figure IA4I shows the fractional 
difference [(Atcode - Atmathematica)/Atmathematica] X 100% for three values of nt and two values of the 
stellar compactness. Figure IA5I shows the fractional difference for three values of eo and ei and two 
values of the stellar compactness. The left and right panels of Figure I A6I show, respectively, the relation 
At(Q;) and the fractional difference given by our code using nt = 100 and eq = ei = 0.01, for four values 
of the stellar compactness that span the range considered in this work. Figures IA4HA6I show that the 
values nt = 100 and eo = ei = 0.01 of these resolution parameters that we use in our waveform code 
when fitting synthetic data provide sufficient accuracy for this purpose. 



A. 1.3. Bolometric flux above the center of an emitting spot on a non-rotating star 

We test the accuracy of our waveform code's computation of the gravitational redshift and deflection 
angle and its integration over the area of the hot spot by comparing the bolometric flux that it predicts 
will be seen by an observer positioned on the radial vector from the center of a non-rotating star through 
the center of the spot, as a function of the compactness of the star and the angular radius of the spot, 
with the analytic result for this flux. 

Consider a uniform circular spot of angular radius A0 on the surface of a non-rotating star with 
compactness MjR, emitting radiation with an isotropic beaming pattern. In the Schwarzschild space- 
time, the bolometric flux seen by a distant observer positioned above the center of the spot is 
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Fig. A5. — Fractional difference between the differential light travel time At as a function of the initial 
ray direction a computed using our code and using Mathematica 8, for three different values of eo and 
ei and two values of the stellar compactness (left: R/M = 5.0; right: R/M = 4.0). This figure shows 
that the values eo = ei = 0.01 used in our waveform fitting are adequate. See Section [A.1.21 for further 
details. 
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Fig. A6. — Left: Differential light travel time At as a function of the initial ray direction a computed 
using our code, for four values of the stellar compactness. The differential light travel time is smaller 
for smaller values of R/M, i.e., for more compact stars. Right: Fractional difference between the light 
travel time At as a function of a computed using our code and using Mathematica 8, for the same four 
values of the stellar compactness. The fractional difference is larger for smaller values of R/M, i.e., for 
more compact stars, but is less than 1.5 x 10~^ for all these values of R/M. This figure shows that the 
resolution parameter values nt = 100 and eg = ei = 0.01 used in our waveform fitting are adequate. We 
have assumed M = 1.4Mq to obtain R from R/M. See Section [A.1.21 for further details. 



where Iq is the specific intensity measured at the stellar surface, \ + z = (1 — 2M/R) ' is the redshift, 
R is the radius of the star, d is the distance to the star, a is the angle that a ray initially makes with 
the surface normal, as measured in the static frame, and ip = ip (a, M/R) is the deflection angle. In 
the limit of vanishing M/R, z — )• 0, cos a -^ cos 9, d cos a/d cos ■0 — >• 1, and we recover the familiar 
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expression 

- / d0 / cos9 dcos9 = TTlo(-) sin^ A0 . (All) 
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A relationship between the angle 9 of the emission point, measured from the symmetry axis of the 
spot, and the angle a a ray makes with the local surface normal is established by the requirement that 
such a ray must be deflected by ■0 = ^ in order to be seen by the observer. In other words, 
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where a{A9) is the value of a for which ip{a) = A9. 

We perform this test of our code as follows. For a given value of M/R and a range of values of 
a, we compute the value of V' that corresponds to each value of a. Next, for each V(q^) (as long as 
V' < vr), we use the relation A0 = V(o^) fo construct the relation a{A9) and then use this value of a 
in expression (IA13P to determine the bolometric flux that is predicted by this analytical expression for 
each value of A0. On the other hand, we use the angular radius A9{a) to compute the bolometric flux 
for each value of a, using our waveform code. We then compare these two values of the flux at each 
value of A0. Figure 1X71 shows an example of such a comparison. The results given by the code differ 
from the results predicted by the analytical expression by less than 5 x 10~^ for A^ > 5°. Even for 
/S.9 < 5°, the difference is less than 7 x 10^'^. 



A. 1.4- Radiation pattern produced by a small spot on a non-rotating star in flat spacetime 

We test the accuracy with which the waveform code computes the radiation pattern from the 
star by using it to compute the flux seen by a distant observer moving around a non-rotating star in 
flat spacetime that has a very small spot that emits isotropically. We then compare this result with 
the analytic result for this situation. In addition to testing the computation of the radiation pattern, 
this comparison also tests the accuracy of the fast Fourier transform (FFT) that we use to represent 
waveforms. 

The radiation pattern for a star and hot spot with the assumed properties is 

/((/>) = maxjcos </>, 0} , (A14) 

where cj) is the phase of the waveform (equivalent to the azimuth of the distant observer, for a non- 
rotating star). The Fourier representation of this waveform is 

oo 

/((/>) = Co + ^c„ cos [n ((/)-(/)„)] , (A15) 

n=l 
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Fig. A7. — Scaled bolometric flux (left) and fractional error in the bolometric flux (right) seen by 
a distant observer, as a function of A0. The agreement between the results given by our waveform 
code (red dots) and the analytic results (blue curve) is excellent. These results are for RjM = 4.0, 
M = IAMq, and 6'spot = ^obs = 0. See Section [A. 1.31 for further details. 
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n is odd (n > 1) 
n is even 



(A16) 



and 



indeterminate (since c„, = 0), n is odd 

0, n = 2,6, 10, , 

vr/n, n = 4, 8, 12, , 



(A17) 



Figures IA8I and IA9I compare these analytic predictions for the fractional rms amplitudes and the 
phases of the first 20 harmonics with the results for these quantities given by our waveform code for 
R/M = 10^, M = 1.4M0, 6'spot = 60°, A(9 = 1°, and e^obs = 20° and the values of the resolution 
parameters listed in Table lAll These figures show that the results given by our waveform code are in 
good agreement with the analytic predictions. 



A. 1.5. Waveform produced by a small spot on a rotating star in flat spacetime 

We test the accuracy with which our waveform code computes the Doppler boost and the photon 
travel time by using it to compute the bolometric waveform produced by a small spot on a rotating star 
in flat spacetime and comparing it with the analytic waveform for this case. 
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Fig. A8. — Fractional rms amplitudes of the first 20 harmonics of the test waveform, showing good 
agreement between the results given by our waveform code and the analytic results. Left: Amplitudes 
given by our waveform code (green dots) compared with the analytic amplitudes (red curve). Right: 
Fractional differences between the amplitudes given by our waveform code and the analytic ampli- 
tudes. Only the differences for harmonic components with nonzero analytic amplitudes are shown. See 
Section [A. 1.41 for further details. 
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Fig. A9. — Phases of the first 20 harmonics of the test waveform, showing good agreement between the 
results given by our waveform code and the analytic results. Left: Phase of the n harmonic component 
of the waveform as a function of harmonic number n (green dots) compared with the analytic phases (red 
curve). Right: Fractional differences between the phases given by our waveform code and the analytic 
phases. Only the results for harmonic components with nonzero analytic amplitudes are shown. See 
Section [A. 1.41 for further details. 



The analytic expression for this waveform is 

f{(p) ex max{5 cosa,0} , 
where (j) is the observed phase (i.e., after correction for the time delay) and 

r/ i\ ^ 

7[1 + (v/c) sin 6'spot sin 6'obs sin $] 



(A18) 



(A19) 
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Fig. AlO. — Top left: Bolometric flux waveform produced by a small spot on the surface of a rotating 
star in flat spacetime, computed using our waveform code (green dots) and the corresponding analytic 
result (red curve), with the time delay set to zero. Top right: Fractional differences (green curve) 
between the fluxes given by our code and the analytic fluxes. Bottom panels: Same as in the top 
panels, but with the time delay included. These results are for R/M = 10^, M/M0=1.6, ^gpot = 60°, 
6»obs = 20°, A6' = 1°, 7^=0.0141357027464 Hz. See Section EXS] for further details. 



is the Doppler factor. Here v is the linear velocity at the rotational equator and ^ = cj) — riAt is the 
phase prior to correction for time delay in terms of the angular rotation frequency i7 of the hot spot 
and the time delay 

(A20) 



At = — (1 — cos a) , 
c 

where a is the angle that the ray makes with the local normal to the stellar surface, as measured in the 

static frame, and is related to other angles via 



cos a = cos 6'spot cos 6'obs + sin e^pot sin 6'obs cos $ . 



(A21) 



To compute the flux at each phase (j) analytically, one must simultaneously solve for <I> and cos a using 
this set of equations. 



Figure lAlOl shows that the bolometric waveform given by our waveform code for this case is in good 
agreement with the analytic result, both when the time delay set to zero and when it is included. 
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A. 1.6. Angular deflection of a pencil beam from a sm.all spot on a rotating star 

We test the accuracy with which our waveform code computes the relativistic aberration and 
gravitational hght-bending effects by computing the rotational phase at which a beam of radiation 
emitted from the stellar surface in a given direction would be seen by a distant observer and comparing 
this result with the phase given by Mathematica 8 (the "predicted" phase). If the star is not rotating, 
the rotational phase is undefined. In this case what is computed is the azimuth at which the beam of 
radiation would be seen by a distant observer. 

In these tests, we assume that the radiation is emitted in a narrow "pencil" beam from a very 
small (A^ = 1°) spot located on the star's rotational equator (^spot = 90°) and that the observer is in 
the plane defined by the star's rotational equator (6'obs = 90°). Because of the spherical symmetry of 
the Schwarzschild spacetime, the radiation will remain in the equatorial plane as it propagates. The 
phase (azimuth) at which the radiation is seen by a distant observer depends on the angle Og of the 
pencil beam relative to the local normal to the stellar surface, as measured in the corotating frame, the 
linear velocity v/c of the hot spot produced by its rotation, which determines the amount of relativistic 
aberration, and the compactness M/R of the star, which determines the amount of light bending. 

We first use a Lorentz transformation to determine the angle oq of the ray in the static frame that 
corresponds to the angle ckq of the ray in the comoving frame. We then compute the angular deflection 
of this ray using Mathematica, which allows us to determine the predicted phase (azimuth) at which 
the beam would be seen by a distant observer. Next we use our waveform code to compute the phase 
(azimuth) at which such a beam would be seen by a distant observer. To do this we set the beaming 
function g{a') equal to 1 for \a' — o'qI < a^i and equal to for |a' — QqI > a^t' , with a^t' = 0.01. 
In order to test the computation of the aberration and the gravitational bending separately from the 
computation of the time delay, we set the light travel time to be zero for all rays. The phase (azimuth) 
at which the beam would be seen by a distant observer is then equal to the total angular deflection of 
the ray produced by the relativistic aberration and light bending. For simplicity, we trace only rays 
that propagate in the prograde direction relative to the rotation of the hot spot. Finally, we compare 
the phase (azimuth) given by the waveform code to the phase predicted using Mathematica. 

Tables IA21 IA31 and IA4I show the results for the computed and predicted phases (or azimuths) at 
which the pencil beam would be seen by a distant observer. These results are for emission from the 
surface of a rapidly rotating star in flat spacetime, from the surface of a non-rotating relativistic star, 
and from the surface of a rapidly rotating relativistic star, for a range of stellar compactnesses and hot 
spot rotation rates. All these results are for M = I.AMq. 



A. 1.7. Observed temperature of thermal emission from a rotating star in flat spacetime, as a function 

of rotational phase 

We test the accuracy with which our waveform code computes the relativistic Doppler effect by 
computing the observed temperature of thermal emission from the surface of a rotating star in flat 
spacetime as a function of the rotational phase and then comparing the result with the analytic result 
for the observed temperature as a function of phase. 

For simplicity, we consider isotropic emission from a small spot at an inclination ^spot on the stellar 
surface, with a Planck spectrum having a temperature Tq as measured in the comoving frame. Then 
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Fig. All. — Left: Temperature of the Planck spectrum seen by a distant observer as a function of 
the star's rotational phase computed using our waveform code (red dots) and the analytic expression 
(blue curve). Right: Fractional difference (green curve) between the observed temperature computed 
using our waveform code and the analytic expression. The results given by our waveform code and the 
analytic results are in excellent agreement. These results are for a star with a rotation rate such that 
the linear velocity at the rotational equator is v/c = 0.7, R/M = 10^, for M = 1.4Mq, ^spot = 60°, 
6(obs = 20°, and Ae* = 1°. See Section [AT3 for further details. 



the spectrum seen by a distant observer with an inclination ^obs is a Planck spectrum with temperature 

To 



n4>) 



(A22) 



7[1 - {vq/c) sin 6lspot sin 6'obs sin </>] 

Also for simplicity, we set the light travel time to zero. The observed phase is then the azimuthal 
position of the spot relative to the azimuth of the observer. To test a range of Doppler boosts, we 
choose values of the stellar rotational frequency that produce a linear velocity at the rotational equator 
of 0.1c, 0.4c, and 0.7c. Figure lA 1 1 1 shows the comparison of the temperature between the code and the 
analytic result, for the case of v/c = 0.7. It is clear that the agreement is excellent. 



A. 2. Tests of the Waveform Fitting Code 

In this subsection, we outline the tests we use to validate the set of codes that fit model waveforms 
to waveform data. An important part of this code set is the code that computes the likelihood of an 
observed waveform, given a set of model parameters. The accuracy of this step depends crucially on the 
accuracy with which the model waveforms are computed. We have shown in the previous subsection 
that our code computes model waveforms with high accuracy. Here we discuss three relatively simple 
end-to-end tests that we use to validate the waveform fitting process. 



A. 2.1. Determining the compactness, temperature, and distance given uniform emission from a 

non-rotating star with no background 

In this exercise, we test our model-fitting code by fitting a model of emission from a non-rotating 
star to the synthetic spectrum produced by a non-rotating star that emits radiation uniformly over its 
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entire surface, with a Planck spectrum. We assume there is no background. The observed spectrum is 
therefore a Planck spectrum with temperature 

Tobs = Tco (l - ^] , (A23) 

where Tco is the temperature of the emission measured in the comoving frame at the stellar surface. 
In this case, the goodness of the fit is determined entirely by the shape and normalization of the 
model spectrum. If the distance is not independently known, the normalization of the spectrum can 
be adjusted by changing the distance. Hence the likelihood of the data given the model has the same 
value (its maximum value) for all Tco-M/R pairs related by equation ()A23p . Consequently, any change 
in the surface temperature of the emission in the model being fit can be perfectly compensated by 
appropriate changes in the compactness M/R of the model star and its distance. This is one of the 
strongest degeneracies in the waveform fitting problem. 

This strong degeneracy can be broken for the synthetic data considered in this example if the mass 
M of the star and the distance d to the star are both known independently of the fitting process, because 
knowledge of these two quantities fixes the normalization of the observed spectrum. To see this, note 
that the normalization of the spectrum is proportional to A/d'^^ where A is the projected area of the 
emitting region. In the case considered here, A = -kE? and is related to the stellar compactness and the 
mass through R = (R/M) x M. Hence if M and d are both known independently, the temperature and 
flux of the observed radiation determines the flux at the stellar surface. This then fixes R and hence 
M/R, since M is assumed known. This result is illustrated in the left-hand panel of Figure IA12I If, 
however, d is not known, the normalization of the spectrum from the stellar surface is unknown and 
can be freely adjusted, even if M is known. In this case, we expect that all models having Tq and 
M/R values that fall on the curve given by equation (IA23P will give good fits. The right-hand panel of 
Figure IA12I illustrates this. 



A. 2. 2. Determining the spot inclination, observer inclination, and distance, given thermal emission 
from a small spot on a non-rotating star in flat spacetime with no background 

In this exercise, we test our model fitting code by fitting a model of emission from a non-rotating 
star in flat spacetime to the synthetic flux produced by a non-rotating star that emits radiation from a 
small spot, with a Planck spectrum. We assume that all the parameters in the model are independently 
known except the inclination of the spot and the observer, or the inclination of the spot and the observer 
and the distance to the star. 

In flat spacetime, if 6'spot + ^obs < f (so that the spot is not occulted by the star), the flux seen by 
an observer at the stellar azimuth (p is 

/((/)) = cos 6'spot cos 6'obs + sin 6'spot sin 6'obs cos (/> . (A24) 

This radiation pattern is a purely sinusoidal function of (j) and can therefore be uniquely specified by 
the mean flux and the fractional rms amplitude of the flux variation with 0. The mean flux is 

/mean = COS 6'spot COS 6'obs , (A25) 

while the fractional rms amplitude of the flux variation with (j) is 

/rms = (I/V2) tan 6'spot tan 6'obs • (A26) 
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Fig. A12. — Code test illustrating the strong degeneracy between the temperature T^o of the emission in 
the comoving frame and the stellar compactness M/R. Left: All parameters in the model except T^o and 
M/R were fixed in advance at the values used to generate the synthetic observed spectrum. The plotted 
points show the values of T^o and R/M that give log likelihoods within 1.15 (i.e., la for 2 degrees of 
freedom) of the log likelihood given by the best-fit model. The models that give high likelihoods cluster 
tightly around T^o = 2.0 keV and R/M = 5.0, the values used in generating the synthetic observed 
spectrum. Right: All parameters except Tco, R/M, and d were fixed in advance at the values used to 
generate the synthetic spectrum. The plotted points show the values of the parameters Tco and R/M 
that give log likelihoods within 1.73 (i.e., la for 3 degrees of freedom) of the log likelihood given by 
the best-fit model. In this case, the models that give high likelihoods cluster tightly along the curve 
given by equation ()A23P with Tobs equal to the observed value, since the normalization of each model 
can be freely adjusted by varying the distance to the star. The results shown here are for M = 1.6Mq, 
/S.9 = 180°, d = 35 kpc, ~ 10^ counts from hot spot, and no background. See Section [A.2.11 for further 
details. 



For a given fractional rms amplitude, equation ()A26P defines a curve in the ^spot-^obs plane along which 
the amplitude remains constant. 

If the distance to the star is independently known, the model will provide a good fit to the data 
only close to two points on the curve given by equation ()A26p . One of these is the point where 0spot 
and ^obs have the values that were used to generate the synthetic observed radiation pattern. The other 
is the point where the values of these two inclinations are interchanged, because the values of /mean 
and /rms are not changed by interchanging the two inclinations. Other points on the curve are strongly 
disfavored, because the value /mean at these points will be very different from the observed value. If 
the distance is not independently known and must be included as a parameter in the fit, good fits can 
be found all along the curve, because the mean observed flux can be adjusted by varying the distance. 
These results are illustrated by the numerical results shown in Figure IA13I 
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Fig. A13. — Code test illustrating the strong degeneracy between the spot and observer inclinations 
^spot and ^obs- Left: All parameters in the model except ^spot and ^obs were fixed in advance at the 
values used to generate the synthetic radiation pattern. The plotted points show the values of ^spot 
and ^obs that give log likelihoods within 1.15 (i.e., la for 2 degrees of freedom) of the log likelihood 
given by the best-fit model. The models that give high likelihoods cluster tightly around two points, 
one corresponding to the values of ^spot = 20° and ^obs = 60° used in generating the synthetic observed 
radiation pattern and the other given by interchanging the two values of ^spot and ^obs) which leaves 
the mean flux, /meam and the rms variation of the flux with cp, frms, unchanged [see equations ()A25P 
and ()A26p ]. Right: All parameters in the model except ^spot and ^obs and d were fixed in advance at the 
values used to generate the synthetic radiation pattern. The plotted points show the values of 6'spot and 
^obs that give log likelihoods within 1.73 (i.e., la for 3 degrees of freedom) of the log likelihood given 
by the best-fit model. In this case, the models that give high likelihoods cluster tightly along the curve 
given by equation (IA26p . where /mean and /rms remain unchanged, but are spread along the curve, 
because the normalization of each model can be freely adjusted by varying the distance to the star. 
The results shown here are for R/M = 10^^ M = 1.6Mq, (9spot = 60°, 6'obs = 20°, A9 = 0.01 radians, 
Tco = 2.0 keV, d = 6 kpc, ~ 10^ counts from the hot spot, and no background. See Section [A. 2. 21 for 
further details. 

A. 2. 3. Determining the compactness, temperature, and distance, given thermal emission from a small 

spot on a non-rotating star in flat spacetime with a background 



This test is the same as the test discussed in Section IA.2.11 except that here we add a background 
that is uniform in azimuth and has a thermal spectrum. In these fits, we treat the background in one 
of two ways: 



1. we add to each model spectrum exactly the background that was used in generating the synthetic 
observed spectrum (i.e., we do not determine the background as part of the fitting procedure), or 

2. we determine the background as part of the fitting procedure, using the code module that imple- 
ments the procedure described in Section 13.3.31 



As subcases, we either set the distance to the value used in generating the synthetic observed spectrum 
during the fitting process (i.e., we assume that the distance is independently known), or we include the 
distance as one of the parameters to be determined in the fitting process. As explained before, assuming 
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Fig. A14. — Code test illustrating the increased uncertainty in the fitted parameters induced by the 
additional fluctuations in the number of counts produced by a background component in the observed 
flux. Left: All parameters in the model except Tco and M/R were fixed in advance at the values used 
to generate the synthetic observed spectrum. The plotted points show the values of Tco and R/M that 
give log likelihoods within 1.15 (i.e., la for 2 degrees of freedom) of the log likelihood given by the 
best-fit model. The models that give high likelihoods cluster around Tco = 2.0 keV and R/M = 5.0, 
the values used in generating the synthetic observed spectrum, but not nearly as tightly as when there 
is no background. Right: All parameters except Tco, R/M, and d were fixed in advance at the values 
used to generate the synthetic spectrum. The plotted points show the values of the parameters Tco and 
R/M that give log likelihoods within 1.73 (i.e., la for 3 degrees of freedom) of the log likelihood given 
by the best-fit model. In this case, the models that give high likelihoods cluster along the curve given 
by equation ()A23P with Tobs equal to the observed value, since the normalization of each model can be 
freely adjusted by varying the distance to the star. Again, they are not nearly as tightly clustered as 
when there is no background. The results shown here are for M = 1.6Mq, IS.6 = 180°, d = 35 kpc, 
~ 10^ counts from the hot spot, Tbkg = 1-5 keV, and ~ 9 x 10^ counts from the background. See 
Section IA.2.31 for further details. 



that the distance is independently known amounts to fixing the normalization of the spectrum, which 
more strongly constrains the fit. 



Figure IA14I shows that the increase in the total number of counts due to the additional counts 
contributed by the background, increases the uncertainties in Tco and M/R. In producing the fits shown 
in this figure, we added to each model spectrum exactly the background that was used in generating 
the synthetic observed spectrum. The left-hand panel in the figure shows the fit assuming that the 
distance to the star is independently known, whereas the right-hand panel shows the fit obtained if the 
distance is determined as part of the fitting procedure. Compared to the case without any background 
(see Figure lA12p . the scatter around the expected relationship between Tco and M/R is larger when a 
background is present that when there is no background, even though the exact background has been 
added to each model. The reason for the increased scatter is the increase in the fluctuation in the total 
number of counts, due to the counts contributed by the background (the fluctuation scales as the square 
root of the total number of counts). 
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Fig. A15. — Percentage differences between three bolometric waveforms computed using three different 
values of the spot resolution parameters Nlat and Nlong and a reference waveform computed using 
Nlat = Nlong = 100, the value we use in our waveform analyses and the code tests described in this 
appendix. This figure shows that for the values of these resolution parameters that we use, the fractional 
error in the bolometric waveform is < 3 x 10~^. See Section [A.3I for further details. 



A. 3. Convergence tests 

We wish to use values for the waveform code resolution parameters listed in Table lAll that will 
provide a resolution fine enough to meet our accuracy requirements, but no finer, so that our code runs 
as fast as possible. In order to determine these values of the resolution parameters, we perform a series 
of convergence tests using different values of these parameters. In these tests, we compare bolometric 
waveforms computed for different values of the resolution parameters, assuming the following values 
of the physical parameters: R/M = 5.0, M/Mq = 1.6, i^rot = 600 Hz, 6*8501 = 90°, A9 = 25°, and 
^obs = 20° . We use the values of the resolution parameters listed in Table lAH except that for clarity in 
the figures we plot the flux at only 16 values of the waveform phase. 



Figure [A15I shows that our choice of 100 grid points in latitude and longitude (Nlat = Nlong = 100) 
resolves the hot spot well enough that the bolometric flux as a function of phase has a fractional error 
< 3 X 10~ , which is adequate for our purposes. Figure [A16I shows that 1000 grid points in the angle 
a' between the ray and the local normal (Nalpha = 1000) resolves this angle well enough that the 
bolometric flux has a fractional error < 10~^, adequate for our purposes. Figure [A17I shows that 1001 
grid points in the final value of the waveform phase (Nloc = 1001) resolves the phase well enough that 
the bolometric flux has a fractional error < 10~^, which is more than adequate for our purposes. Indeed, 
Nloc = 101 provides essentially the same accuracy. 

The tests of the values of the resolution parameters used in computing the deflection angle and the 
light travel time have been presented in Figures IA1IIA21 \A4\ and IA5[ 
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Fig. A16. — Percentage differences between tliree bolometric waveforms computed using three different 
values of the ray direction resolution parameter Nalpha and a reference waveform computed using 
Nalpha = 1000, the value we use in our waveform analyses and the code tests described in this appendix. 
This figure shows that for the values of these resolution parameters that we use, the fractional error in 
the bolometric waveform is < 10~^. See Section [A. 31 for further details. 
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Fig. A17. — Percentage differences between three bolometric waveforms computed using three different 
values of Nloc, the grid spacing in the waveform phase with time delays included, and a reference 
waveform computed using Nloc = 1001, the value we use in the code tests described in this appendix. 
This figure shows that the fractional error in the bolometric waveform for this value of Nloc is < 10~^. 
See Section [A. 31 for further details. 
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Table Al. Values of the waveform code resolution parameters 

Parameter Value 

Nlat (spot grid points in stellar latitude) 100 

Nlong (spot grid points in stellar longitude) 100 

Nalpha (grid points in a) 1000 

Nphase (grid points in phase w. time delays included) 1000 

Nloc (grid points in phase w.o. time delays included) 1001 

Ud (grid points used in computing the deflection angle) 100 

ed (deflection angle integration parameter) 0.1 

nt (grid points used in computing the time delay) 100 

€o (first time delay integration parameter) 0.01 

€i (second time delay integration parameter) 0.01 
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Table A2. Comparison of pencil-beam phases'^ 



ao 


v/c 


phase (code) 


phase (predicted) 


15° 


0.1 


0.958 


0.958 


30° 


0.1 


0.917 


0.917 




0.4 


0.917 


0.917 


45° 


0.1 


0.875 


0.875 




0.4 


0.875 


0.875 




0.7 


0.877 


0.875 


60° 


0.1 


0.833 


0.833 




0.4 


0.833 


0.833 




0.7 


0.833 


0.833 


75° 


0.1 


0.792 


0.792 




0.4 


0.791 


0.792 




0.7 


0.791 


0.792 


90° 


0.1 


0.751 


0.750 




0.4 


0.751 


0.750 




0.7 


0.751 


0.750 



'^Comparison of the computed and predicted 
phases at which a narrow beam of radiation emit- 
ted from a small spot on the equator of a rotating 
star at an angle ao relative to the local surface 
normal (as measured in the static frame) is seen 
by a distant observer who is also located in the 
equatorial plane. In this test, R/M = 10^, so 
the spacetime is essentially flat. The stellar rota- 
tion frequency is chosen so the linear velocity of 
the stellar surface in the rotation equator is 0.1c, 
0.4c, or 0.7c. The values of v/c that have been 
omitted from the table would produce beams that 
point in the retrograde direction in the comoving 
frame, a geometry not considered in this test. See 
Section [A. 1 .61 for further details. 
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Table A3. Comparison of pencil-beam phases'' 



ao 


R/M 


azimuth (code) 


azimuth (predicted) 


0° 


4.0 


0.998 


1.000 




5.0 


0.998 


1.000 




6.0 


0.999 


1.000 


15° 


4.0 


0.941 


0.941 




5.0 


0.946 


0.946 




6.0 


0.948 


0.949 


30° 


4.0 


0.881 


0.881 




5.0 


0.891 


0.892 




6.0 


0.897 


0.897 


45° 


4.0 


0.818 


0.818 




5.0 


0.835 


0.836 




6.0 


0.844 


0.845 


60° 


4.0 


0.751 


0.751 




5.0 


0.777 


0.777 




6.0 


0.790 


0.790 


75° 


4.0 


0.674 


0.674 




5.0 


0.713 


0.713 




6.0 


0.733 


0.733 


90° 


4.0 


0.579 


0.576 




5.0 


0.642 


0.640 




6.0 


0.671 


0.669 



'^Comparison of the computed and predicted azimuths 
at which a narrow beam of radiation emitted from a 
small spot on the equator of a non-rotating star at an 
angle oq relative to the local surface normal (as measured 
in the static frame) would be seen by a distant observer 
who is also located in the equatorial plane, for three 
values of the stellar compactness. See Section IA.1.61 for 
details. 
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Table A4. Comparison of pencil-beam phases'' 



ao 


R/M 


v/c 


phase (code) 


phase (predicted) 


30° 


3.5 


0.4 


0.870 


0.871 




6.0 


0.4 


0.897 


0.897 


45° 


3.5 


0.4 


0.801 


0.802 






0.7 


0.802 


0.802 




6.0 


0.4 


0.845 


0.845 






0.7 


0.790 


0.790 


60° 


3.5 


0.4 


0.725 


0.725 






0.7 


0.725 


0.725 




6.0 


0.4 


0.790 


0.790 






0.7 


0.790 


0.790 


75° 


3.5 


0.4 


0.632 


0.633 






0.7 


0.632 


0.633 




6.0 


0.4 


0.733 


0.733 






0.7 


0.732 


0.733 


90° 


3.5 


0.4 


0.498 


0.495 






0.7 


0.498 


0.495 




6.0 


0.4 


0.671 


0.669 






0.7 


0.670 


0.669 



'^Comparison of the computed and predicted phases at 
which a narrow beam of radiation emitted from a small 
spot on the equator of a rotating star at an angle uq rela- 
tive to the local surface normal (as measured in the static 
frame) is seen by a distant observer who is also located 
in the equatorial plane, for two values of the stellar com- 
pactness. The stellar rotation frequency is chosen so the 
linear velocity of the stellar surface in the rotation equa- 
tor is either 0.4c or 0.7c. The values of v/c that have been 
omitted from the table would produce beams that point 
in the retrograde direction in the comoving frame, a ge- 
ometry not considered in this test. See Section IA.1.61 for 
further details. 
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B. JOINT FITS 

In this Appendix we discuss the constraints on the waveform model parameters that can be obtained 
by jointly fitting a collection of different, independent sets of waveform data. These sets could be 
different segments of a single burst or different bursts from the same star. 

We show here that the constraints on the interesting parameters that can be obtained by jointly 
fitting many data sets are often comparable to the constraints that would have been obtained by fitting 
a single data set that has an oscillation profile which is the same as the average oscillation profile of 
the multiple data sets and has the same total number of counts as the multiple data sets have when 
combined. This is true even though one might intuitively expect that the extra parameters required to 
fit many data sets in which the uninteresting parameters change from set to set would compromise the 
constraints on the (unchanging) interesting parameters. 

We assume that the values of some of the parameters in the model waveform, such as the mass 
M and radius R of the star, the inclination ^obs of the observer, and the distance d to the star remain 
unchanged during a burst and from burst-to-burst, whereas the values of other parameters, such as the 
angular radius A9 of the hot spot, the inclination 6c of the spot center, and the color temperature Tco 
of the emission from the spot, may change during a burst and from burst-to-burst. The parameters 
that remain unchanged must be treated differently from the parameters that may change. 

Suppose that we have several independent data sets and wish to extract joint constraints on the 
subset of waveform parameters «» that have the same values in all the data sets. We denote the other 
waveform parameters, which may change between data sets, by 13 j. We can derive the correct approach 
to jointly fitting multiple data sets by starting from the Bayesian expression for the unnormalized 
posterior probability density 

qia^) = I dPjP{a,,pj)C{a,,pj) . (Bl) 

Here P is the joint prior probability density over all the parameters and C is the likelihood of all the 
data, given the model and specific values for all the model parameters. If there are multiple data sets, 
then the joint likelihood is the product of the individual likelihoods, i.e., 

£ = []£fc(a„/3,-fc), (B2) 

k 

where here we write Pj^ to indicate the set of potentially variable parameters that is being evaluated 
for data set k. 

We assume that the prior probability density distribution for the fixed parameters «» is the same 
for each data set. Then we can write the full prior probability distribution as 

P{ai,pjl,pj2,■■■)=P{a^)Pl{PJl\a^)P2{PJ2\a^)■■■ (B3) 

where the product is over all data sets and the vertical bar indicates a conditional probability. For 
example, Pi(/3ji|aj) is the prior probability of the variable parameters /3ji in the first data set given 
that the fixed parameters have the values a,. 

The posterior probability density of our fixed parameters is thus proportional to 



q{ai) oc P(ai)]^ 



d(3jkPk{(3jk\ai)£.k{ai, (3jk) 



(B4) 
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Thus the posterior probabihty density is equal to the prior probabihty density times the product of 
the integrals over the variable parameters for each data set, of the prior probability density of the 
variable parameters times the likelihood of the data given the values of all the parameters. The process 
of integrating over the uninteresting parameters (called nuisance parameters in this context) is called 
marginalization. Marginalization over the nuisance parameters means that in some sense they do 
not count in the tally of parameters that are adjusted during the fitting process. It is this part of 
the procedure that ultimately yields constraints from many data sets that are not muddled by the 
introduction of many fit parameters. 

It is important to note that the posterior probability density is usually not the product of the 
marginalized probability densities Qk{cti) for each of the data sets. That is, 



/(«.) / n 



dl3jkPiai)Pil3jk\ai)Ckiai, /3jk) 



(B5) 



k 

For n data sets, this incorrect expression for the posterior probability density would yield 



g(a,)oc[P(a,)rn 



dl3jkPk{l3jk)Ckiai, I3jk) 



(B6) 



In this incorrect expression for n data sets, there are n factors of the prior probability density, compared 
with only one in the correct expression. To see that n factors is incorrect, imagine a case in which the 
prior probability distribution is a multivariate Gaussian and that only one of a large number n of data 
sets is informative. The incorrect expression would nevertheless produce very tight constraints, because 
it raises the Gaussian to a high power. 



The posterior probability density ()B4p is not in general proportional to the posterior probability 
density one would obtain from an equivalent number of counts in a profile that has the average shape of 
the profile in the multiple data sets. To see this, let's perform a thought experiment in which we have 
obtained a single, continuous segment of data during which all the waveform parameters are known to 
be fixed, but we analyze it as some number of contiguous data sets where the values of the parameters 
Oi are fixed but the values of the Pjk need not be. Then the joint analysis gives equation (IB4p . whereas 
the analysis that assumes all the parameters are fixed gives 

q{ai) ex P{ai) f dl3jP{(3j\ai)C(ai, (3j) , (B7) 

where as before C = Y\kJ~-k{oii, (3j). Note that instead of allowing the /3j to take on different values 13 jk 
in each data set, here they are assumed to have the same values. In general, equation ()B7p does not 
give the same posterior as equation ()B4p . 



There are, however, circumstances in which the posteriors are identical. Consider a situation in 

which (1) the prior for the variable parameters is independent of the fixed parameters, and (2) the 
likelihood is the product of two independent functions, one for ai and one for fSjk'- 

^k{oii,/3jk) = Cka^kl3 ■ (B8) 

Then 

q{ai) oc P{ai) Hfc [/ dPjkPk{l3jk)^k{ai, /3jk)] .g„x 

= [^(«0 Uk ^ka] Uk U dP,kPk{/3jk)Ckp] . ^ ' 
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Fig. Bl. — An example showing that the constraints on M and R obtained by jointly analyzing multiple 
data sets as described in the text are usually comparable to the constraints obtained by analyzing a 
single data set with an oscillation profile similar to the average profile of the oscillations in the multiple 
data sets and with the same total number of counts as there are in the multiple data sets. The results 
shown here are for the high-inclination reference case with our medium background and no independent 
knowledge of the parameters in the model. Left: Constraints obtained by analyzing a single data set 
with about 10^ counts from the hot spot. Right: Constraints obtained by breaking up the single data 
set into five segments of equal length and then analyzing the segments jointly. 



The second factor in equation (IBOh is independent of a and therefore enters the posterior as a constant, 
which will be normalized away. The product W^ Cka is just the total likelihood, and the posterior will 
therefore be unchanged by dividing the original single data set into multiple data sets. 

In order to investigate whether this general result applies to the problem studied in this paper, 
namely, analysis of multiple segments of burst oscillation data, we have carried out a number tests. In 
these tests we compared the constraints obtained by analyzing a set of waveform data as single segment 
and as a sequence of segments. We find that in practice the constraints obtained by jointly analyzing 
multiple data sets are usually comparable to the constraints obtained by analyzing a single data set 
with an oscillation profile similar to the average profile of the oscillations in the multiple data sets and 
with the same total number of counts as there are in the multiple data sets. 

We illustrate these results by the following two examples. Both compare the constraints on M 
and R obtained by analyzing a single data set that has about 10® counts from the hot spot with the 
constraints obtained by breaking up the single data set into five segments of equal length and then 
analyzing the segments jointly, using the method described above. All the confidence regions shown 
here are for the high-inclination reference case and our medium background. Figure IBll compares the 
constraints obtained when no independent information about the parameters in the model is available. 
Figure IB2I compares the constraints obtained when the distance to the star is known independently 
of the waveform analysis. In both cases, the constraints on M and R obtained by jointly analyzing 
multiple data sets as described in the text are comparable to the constraints obtained by analyzing a 
single data set with an oscillation profile similar to the average profile of the oscillations in the multiple 
data sets and with the same total number of counts as there are in the multiple data sets. 
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Fig. B2. — Another example showing that the constraints on M and R obtained by jointly analyzing 
multiple data sets as described in the text are usually comparable to the constraints obtained by 
analyzing a single data set with an oscillation profile similar to the average profile of the oscillations in 
the multiple data sets and with the same total number of counts as there are in the multiple data sets. 
The results shown here are for the high-inclination reference case with our medium background and 
independent knowledge of the distance to the star. Left: Constraints obtained by analyzing a single 
data set with about 10® counts from the hot spot. Right: Constraints obtained by breaking up the 
single data set into five segments of equal length and then analyzing the segments jointly. 
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